sapply(c(
         "rjson", 
         "data.table", 
         "dplyr", 
         "ggplot2", 
         "stringr", 
         "purrr", 
         "foreach", 
         "doParallel", 
         "patchwork", 
         "lme4", 
         "lmerTest",
         "testit",
         "latex2exp",
         "ggpubr"
         ), 
       require, character=TRUE)
## Loading required package: rjson
## Loading required package: data.table
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:data.table':
## 
##     between, first, last
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## Loading required package: ggplot2
## Loading required package: stringr
## Loading required package: purrr
## 
## Attaching package: 'purrr'
## The following object is masked from 'package:data.table':
## 
##     transpose
## Loading required package: foreach
## 
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
## 
##     accumulate, when
## Loading required package: doParallel
## Loading required package: iterators
## Loading required package: parallel
## Loading required package: patchwork
## Loading required package: lme4
## Loading required package: Matrix
## Loading required package: lmerTest
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
## Loading required package: testit
## Loading required package: latex2exp
## Loading required package: ggpubr
##      rjson data.table      dplyr    ggplot2    stringr      purrr    foreach 
##       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE 
## doParallel  patchwork       lme4   lmerTest     testit  latex2exp     ggpubr 
##       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE
sf <- function() sapply(paste0("./Functions/", list.files("./Functions/", recursive=TRUE)), source) # Source all fxs
sf()
##         ./Functions/General/FindDelay.R ./Functions/General/FindDelayAlt.R
## value   ?                               ?                                 
## visible FALSE                           FALSE                             
##         ./Functions/General/Utils.R ./Functions/Model/Models/RunMBayesLearner.R
## value   ?                           ?                                          
## visible FALSE                       FALSE                                      
##         ./Functions/Model/Models/RunMBayesLearnerDiffBeta.R
## value   ?                                                  
## visible FALSE                                              
##         ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayes.R
## value   ?                                                             
## visible FALSE                                                         
##         ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayesAlt.R
## value   ?                                                                
## visible FALSE                                                            
##         ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayes.R
## value   ?                                                          
## visible FALSE                                                      
##         ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayesAlt.R
## value   ?                                                             
## visible FALSE                                                         
##         ./Functions/Model/Models/RunMQLearner.R
## value   ?                                      
## visible FALSE                                  
##         ./Functions/Model/Models/RunMQLearnerDecayTo0Inits.R
## value   ?                                                   
## visible FALSE                                               
##         ./Functions/Model/Models/RunMQLearnerDecayToPessPrior.R
## value   ?                                                      
## visible FALSE                                                  
##         ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffBeta.R
## value   ?                                                              
## visible FALSE                                                          
##         ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffLR.R
## value   ?                                                            
## visible FALSE                                                        
##         ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPrior.R
## value   ?                                                          
## visible FALSE                                                      
##         ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorESAndEpsFixed.R
## value   ?                                                                       
## visible FALSE                                                                   
##         ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorNoCK.R
## value   ?                                                              
## visible FALSE                                                          
##         ./Functions/Model/ModelUtils/AltPlotSimTrainingCurves.R
## value   ?                                                      
## visible FALSE                                                  
##         ./Functions/Model/ModelUtils/aModelUtils.R
## value   ?                                         
## visible FALSE                                     
##         ./Functions/Model/ModelUtils/PlotAllWorstTestOptionsSimVsEmp.R
## value   ?                                                             
## visible FALSE                                                         
##         ./Functions/Model/ModelUtils/PlotSimEmpTest.R
## value   ?                                            
## visible FALSE                                        
##         ./Functions/Model/ModelUtils/PlotSimEmpTestNoSim.R
## value   ?                                                 
## visible FALSE                                             
##         ./Functions/Model/ModelUtils/PrepDfForModel.R
## value   ?                                            
## visible FALSE                                        
##         ./Functions/Model/ModelUtils/PrepDfForModelPROpt.R
## value   ?                                                 
## visible FALSE                                             
##         ./Functions/Model/ModelUtils/RecodeDfs.R
## value   ?                                       
## visible FALSE                                   
##         ./Functions/Model/ModelUtils/SimplePlotSimTrainingDelay.R
## value   ?                                                        
## visible FALSE                                                    
##         ./Functions/Model/OptFxs/RunOptTrainSIT.R
## value   ?                                        
## visible FALSE                                    
##         ./Functions/Model/OptFxs/RunOptTrainSITPR.R
## value   ?                                          
## visible FALSE                                      
##         ./Functions/Model/SimFxs/RunSimTrainSIT.R
## value   ?                                        
## visible FALSE                                    
##         ./Functions/Model/SimFxs/RunSimTrainSITForPR.R
## value   ?                                             
## visible FALSE                                         
##         ./Functions/Model/TrialWiseComps/aModelFunctions.R
## value   ?                                                 
## visible FALSE                                             
##         ./Functions/Model/TrialWiseComps/CalcBayesMean.R
## value   ?                                               
## visible FALSE                                           
##         ./Functions/Model/TrialWiseComps/CalcBayesStd.R
## value   ?                                              
## visible FALSE                                          
##         ./Functions/Model/TrialWiseComps/CalcBayesVar.R
## value   ?                                              
## visible FALSE                                          
##         ./Functions/Model/TrialWiseComps/DecayChoiceKernel.R
## value   ?                                                   
## visible FALSE                                               
##         ./Functions/Model/TrialWiseComps/DecayQVals.R
## value   ?                                            
## visible FALSE                                        
##         ./Functions/Model/TrialWiseComps/UpdateABMatrix.R
## value   ?                                                
## visible FALSE
DefPlotPars()
registerDoParallel(cores=round(detectCores()*2/3))

Load data from studies 1 and 2

s1_train <- read.csv("../data/cleaned_files/s1_train_with_delay_deident.csv")
s1_sit <- read.csv("../data/cleaned_files/s1_SIT_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s2_train <- read.csv("../data/cleaned_files/s2_train_with_delay_deident.csv")
s2_sit <- read.csv("../data/cleaned_files/s2_sit_deident_corrected_names.csv")
s1_pst <- read.csv("../data/cleaned_files/s1_PST_deident.csv") %>% rename(ID=deident_ID)
s2_pst <- read.csv("../data/cleaned_files/s2_PST_deident.csv") %>% rename(ID=deident_ID)

Harmonize variables and create some separate vars

# Study 2 harmonize  
s2_sit[s2_sit$condition=="cogn", "condition"] <- "cognitive" 
s2_train[s2_train$trial_within_condition <= 20, "block"] <- 1
s2_train[s2_train$trial_within_condition > 20, "block"] <- 2

s1_sit$probability <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 1)))
s1_sit$valence <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 2)))

s2_sit$probability <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 1)))
s2_sit$valence <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 2)))

Define paths and functions

# MLE results  
allp <- "../model_res/opts_mle_paper_final/all/"
# Empirical Bayes results 
bp <- "../model_res/opts_mle_paper_final/best/"
# Simulations  
sp <- "../model_res/sims_clean/sims_from_mle/"
# Read model function 
rm <- function(path, model_str) read.csv(paste0(path, model_str))

Plotting delay

s1_train_delay_summs <- s1_train %>% filter(!is.na(delay_trunc)) %>% group_by(condition, delay_trunc, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition', 'delay_trunc'. You can
## override using the `.groups` argument.
s1_tr_delay_summs <- 
  Rmisc::summarySEwithin(s1_train_delay_summs,
                        measurevar = "m",
                        withinvars = c("condition", "delay_trunc"), 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition, delay_trunc
s2_train_delay_summs <- s2_train %>% filter(!is.na(delay_trunc)) %>% group_by(condition, delay_trunc, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition', 'delay_trunc'. You can
## override using the `.groups` argument.
s2_tr_delay_summs <- 
  Rmisc::summarySEwithin(s2_train_delay_summs,
                        measurevar = "m",
                        withinvars = c("condition", "delay_trunc"), 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition, delay_trunc
a <- ggplot(s1_tr_delay_summs, aes(x=delay_trunc, y=m, group=condition, color=condition)) + 
  geom_line(size=2, position = position_dodge(width = .2)) +
  geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3, color="gray57") + 
  geom_hline(yintercept = .5, size=2, color="black") +
  geom_errorbar(aes(x=delay_trunc, y=m, ymin=m-se, ymax=m+se, group=condition),
                position = position_dodge(width = .2), color="black", size=1.5, width=.15) +
  geom_point(size=6, position = position_dodge(width = .2), pch=21, fill="gray57", alpha=.95) +
  ga + ap + ft + tol + ylim(.48, .82) + tp +
  scale_color_manual(values=c("purple", "orange"), 
                     labels=c("cognitive", "overt"))  + ylab("Proportion correct") + xlab("Delay") + 
  ggtitle("Experiment 1")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
b <- ggplot(s2_tr_delay_summs, aes(x=delay_trunc, y=m, group=condition, color=condition)) + 
  geom_line(size=2, position = position_dodge(width = .2)) +
  geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3, color="gray57") + 
  geom_hline(yintercept = .5, size=2, color="black") +
  geom_errorbar(aes(x=delay_trunc, y=m, ymin=m-se, ymax=m+se, group=condition),
                position = position_dodge(width = .2), color="black", size=1.5, width=.15) +
  geom_point(size=6, position = position_dodge(width = .2), pch=21, fill="gray57", alpha=.95) +
  ga + ap + ft + lp + ylim(.48, .82) + tp +
  scale_color_manual(values=c("purple", "orange"), 
                     labels=c("cognitive", "overt"))  + ylab("") + xlab("Delay") + 
  ggtitle("Experiment 2")
delay_comb_plot <- a+b
delay_comb_plot

#ggsave("../paper/figs/fig_parts/delay_plot.png", delay_comb_plot, height=7, width=10, dpi=700)

Supplemental material model validation (although learning curves exp 2 plot is in key_results)

6/5/23 — reran sims to mkae sure completely updated/bug fixed version of par results used

s1_train_sim_m1_eb <- 
  rm(sp, 
     "SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDiffDecayToPessPrior57088.csv")

s2_train_sim_m1_eb <- 
  rm(sp, 
    "SIM_EMPIRICAL_BAYES_study_2_train_SIT__train_RunMQLearnerDiffDecayToPessPrior28414.csv")

s1_test_sim_m1_eb <- 
  rm(sp, 
     "SIM_EMPIRICAL_BAYES_study_1_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior57088.csv")

s2_test_sim_m1_eb <- 
  rm(sp, 
    "SIM_EMPIRICAL_BAYES_study_2_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior28414.csv")

Capture of delay effects

out_a <- SimplePlotSimTrainingDelay(emp_df=s1_train, sim_df=s1_train_sim_m1_eb)
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'iter'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## Warning: Duplicated aesthetics after name standardisation: colour
out_b <- SimplePlotSimTrainingDelay(emp_df=s2_train, sim_df=s2_train_sim_m1_eb) 
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'iter'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## Warning: Duplicated aesthetics after name standardisation: colour
out_a
## Warning: Duplicated aesthetics after name standardisation: colour

out_b
## Warning: Duplicated aesthetics after name standardisation: colour

# ggsave("../paper/figs/fig_parts/simvsemp_delay_plot_1.png", out_a, height=8, width=12, dpi=700)
# ggsave("../paper/figs/fig_parts/simvsemp_delay_plot_2.png", out_b, height=8, width=12, dpi=700)

Worst-option Proportions in Test Phase and comparison to same model without CK

worst_plots_a_ck <- PlotAllWorstTestOptionsSimVsEmp(s1_sit, s1_test_sim_m1_eb, m_string="Simulations \nwith Choice Kernel")
w_a_ck <- worst_plots_a_ck[[1]]+worst_plots_a_ck[[2]]
w_a_ck 

worst_plots_b_ck <- PlotAllWorstTestOptionsSimVsEmp(s2_sit, s2_test_sim_m1_eb, m_string="Simulations \nwith Choice Kernel")
w_b_ck <- worst_plots_b_ck[[1]]+worst_plots_b_ck[[2]]
w_b_ck
## Warning: Removed 1 rows containing missing values (`geom_point()`).

# ggsave("../paper/figs/fig_parts/CK_simvsemp_worst_plot_1.png", w_a_ck, height=6, width=15, dpi=700)
# ggsave("../paper/figs/fig_parts/CK_simvsemp_worst_plot_2.png", w_b_ck, height=6, width=15, dpi=700)

Comparative worst plots

s1_test_sim_m1_eb_no_ck <- 
  rm(sp, 
     "SIM_EMPIRICAL_BAYES_study_2_train_SIT__sit_RunMQLearnerDiffDecayToPessPriorNoCK16278.csv")

s2_test_sim_m1_eb_no_ck <- 
  rm(sp, 
    "SIM_EMPIRICAL_BAYES_study_1_train_SIT__sit_RunMQLearnerDiffDecayToPessPriorNoCK31508.csv")
worst_plots_nck_a <- PlotAllWorstTestOptionsSimVsEmp(s1_sit, s1_test_sim_m1_eb_no_ck, m_string="Simulations \nNo Choice Kernel")
w_a_nck <- worst_plots_nck_a[[1]]+worst_plots_nck_a[[2]]
w_a_nck

worst_plots_nck_b <- PlotAllWorstTestOptionsSimVsEmp(s2_sit, s2_test_sim_m1_eb_no_ck, m_string="Simulations \nNo Choice Kernel")

Captures the mis-spec

w_b_nck <- worst_plots_nck_b[[1]]+worst_plots_nck_b[[2]]
w_b_nck

# ggsave("../paper/figs/fig_parts/NOCK_simvsemp_worst_plot_1.png", w_a_nck, height=6, width=15, dpi=700)
# ggsave("../paper/figs/fig_parts/NOCKsimvsemp_worst_plot_2.png", w_b_nck, height=6, width=15, dpi=700)

Parameter recovery

m1_par_recov <- 
  read.csv("../model_res/par_recov_clean/pr_opts_mle/best_pr/par_recov_BEST_EMPIRICAL_BAYES_mv_gaussstudy_1_train_SIT_RunMQLearnerDiffDecayToPessPrior_76217.csv") # Confirmed created 3/28 after the 3/23 bug fix on decay  
eb_recov_pars_m1 <- m1_par_recov[grep("EB_recovered", names(m1_par_recov))]
simmed_pars_m1 <- m1_par_recov[grep("simmed", names(m1_par_recov))]

Phi difference

this_eb_pr_df <- 
  data.frame("simulated"=simmed_pars_m1$cog_phi_simmed-simmed_pars_m1$overt_phi_simmed,
             "EB_recovered"=eb_recov_pars_m1$cog_phi_EB_recovered-eb_recov_pars_m1$overt_phi_EB_recovered)

phi_diff <- ggplot(this_eb_pr_df, aes(x=simulated, y=EB_recovered)) +
      geom_line(aes(x=simulated, y=simulated), size=4) +
      geom_point(size=8, pch=21, fill="gray57", alpha=.8) + 
      stat_cor(method="spearman", size=6, label.y=.9) +
      ggtitle(TeX("$\\phi^{Cog}-\\phi^{Overt}")) +
      ga + ap + tp  + 
      ylab("Recovered") + xlab("Simulated")
phi_diff

#ggsave("../paper/figs/fig_parts/pr/pr_phi_diff.png", phi_diff, height=4, width=8, dpi=700)

Percent correctly recovering the sign difference

this_eb_pr_df$sim_sign <- if_else(this_eb_pr_df$simulated >= 0, 1, 0)
this_eb_pr_df$eb_sign <- if_else(this_eb_pr_df$EB_recovered >= 0, 1, 0)
# this_eb_pr_df$sim_sign==this_eb_pr_df$eb_sign
# (this_eb_pr_df$sim_sign==this_eb_pr_df$eb_sign)*1
ta <- table((this_eb_pr_df$sim_sign==this_eb_pr_df$eb_sign)*1)
ta[2]/sum(ta)
##         1 
## 0.7142857
plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$epsilon_simmed, eb_recov_pars_m1$epsilon_EB_recovered, 
                     "$\\epsilon", stat_pos = .32)
plot

# ggsave("../paper/figs/fig_parts/pr/eps.png",
#        plot, height=4, width=8, dpi=700)
plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$cog_phi_simmed, eb_recov_pars_m1$cog_phi_EB_recovered, 
                     "$\\phi^{Cog}")
plot

# ggsave("../paper/figs/fig_parts/pr/phi_cog.png",
#        plot, height=4, width=8, dpi=700)
plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$overt_phi_simmed, eb_recov_pars_m1$overt_phi_EB_recovered, 
                     "$\\phi^{Overt}")
plot

# ggsave("../paper/figs/fig_parts/pr/phi_overt.png",
#        plot, height=4, width=8, dpi=700)
plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$choice_LR_simmed, eb_recov_pars_m1$choice_LR_EB_recovered, 
                     "$\\alpha_{CK}")
plot

# ggsave("../paper/figs/fig_parts/pr/ck.png",
#        plot, height=4, width=8, dpi=700)
plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$explor_scalar_simmed, eb_recov_pars_m1$explor_scalar_EB_recovered, 
                     "ES")
plot

# ggsave("../paper/figs/fig_parts/pr/es.png",
#        plot, height=4, width=8, dpi=700)
plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$q_LR_simmed, eb_recov_pars_m1$q_LR_EB_recovered, 
                     "$\\alpha_{Q}")
plot

# ggsave("../paper/figs/fig_parts/pr/qlr.png",
#        plot, height=4, width=8, dpi=700)
plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$beta_simmed, eb_recov_pars_m1$beta_EB_recovered, 
                     "$\\beta", stat_pos = 25)
plot

Simulation with a high difference in phi

Load model fits

m1_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior17352.csv")
m1_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior58319.csv")

m1_study1_eb <- rbind(m1_study1_eb_v1, m1_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m1_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior12123.csv")
m1_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior74336.csv")

m1_study2_eb <- rbind(m1_study2_eb_v1, m1_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

Find a participant with relatively high diff in decay

m1_study1_eb$cog_phi[25]
## [1] 0.6817201
m1_study1_eb$overt_phi[25]
## [1] 0.1300788
m1_study1_eb[25, ]
## # A tibble: 1 × 12
## # Groups:   ID [1]
##       X epsilon  q_LR cog_phi overt_phi choice_LR explor_scalar  beta
##   <int>   <dbl> <dbl>   <dbl>     <dbl>     <dbl>         <dbl> <dbl>
## 1    25  0.0365 0.732   0.682     0.130    0.0425         0.220  4.97
## # ℹ 4 more variables: convergence <int>, nll <dbl>, AIC <dbl>, ID <int>

This did correspond to a big difference in performance empirically

s1_train %>% filter(ID==25) %>% group_by(valence_and_probability, condition) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'valence_and_probability'. You can override
## using the `.groups` argument.
## # A tibble: 8 × 3
## # Groups:   valence_and_probability [4]
##   valence_and_probability condition     m
##   <chr>                   <chr>     <dbl>
## 1 40-10_punishment        cognitive 0.6  
## 2 40-10_punishment        overt     0.95 
## 3 40-10_reward            cognitive 0.9  
## 4 40-10_reward            overt     1    
## 5 90-10_punishment        cognitive 0.525
## 6 90-10_punishment        overt     0.775
## 7 90-10_reward            cognitive 0.05 
## 8 90-10_reward            overt     1
s1_sit %>% filter(ID==25) %>% group_by(valence_and_probability, condition) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'valence_and_probability'. You can override
## using the `.groups` argument.
## # A tibble: 8 × 3
## # Groups:   valence_and_probability [4]
##   valence_and_probability condition     m
##   <chr>                   <chr>     <dbl>
## 1 40-10_punishment        cognitive  0.75
## 2 40-10_punishment        overt      1   
## 3 40-10_reward            cognitive  1   
## 4 40-10_reward            overt      1   
## 5 90-10_punishment        cognitive  1   
## 6 90-10_punishment        overt      1   
## 7 90-10_reward            cognitive  0   
## 8 90-10_reward            overt      1

Summarize the correct and incorrect q-values

corr_q_vals_summ <- 
  s1_train_sim_m1_eb %>% filter(ID==25) %>% 
  group_by(cond, probability, valence, trial_within_condition) %>% summarize(m=mean(sim_correct_Q_vals))
## `summarise()` has grouped output by 'cond', 'probability', 'valence'. You can
## override using the `.groups` argument.
corr_q_vals_summ$valence <- factor(corr_q_vals_summ$valence, levels=c("reward", "punishment"))
corr_q_vals_summ[corr_q_vals_summ$cond=="cog", "cond"] <- "cognitive"
incorr_q_vals_summ <- 
  s1_train_sim_m1_eb %>% filter(ID==25) %>% 
  group_by(cond, probability, valence, trial_within_condition) %>% summarize(m=mean(sim_incorrect_Q_vals))
## `summarise()` has grouped output by 'cond', 'probability', 'valence'. You can
## override using the `.groups` argument.
incorr_q_vals_summ$valence <- factor(incorr_q_vals_summ$valence, levels=c("reward", "punishment"))

incorr_q_vals_summ[incorr_q_vals_summ$cond=="cog", "cond"] <- "cognitive"

Summarize the simulated performance during learning

sim_perf_summ <- 
  s1_train_sim_m1_eb %>% filter(ID==25) %>% 
  group_by(cond, probability, valence, trial_within_condition) %>% summarize(m=mean(sim_corrs))
## `summarise()` has grouped output by 'cond', 'probability', 'valence'. You can
## override using the `.groups` argument.
sim_perf_summ$valence <- factor(sim_perf_summ$valence, levels=c("reward", "punishment"))

sim_perf_summ[sim_perf_summ$cond=="cog", "cond"] <- "cognitive"

Take the difference in q-values

assert(all(corr_q_vals_summ$trial_within_condition==incorr_q_vals_summ$trial_within_condition))
assert(all(corr_q_vals_summ$valence==incorr_q_vals_summ$valence))
assert(all(corr_q_vals_summ$probability==incorr_q_vals_summ$probability))

corr_q_vals_summ$diff <- corr_q_vals_summ$m-incorr_q_vals_summ$m
a <- 
  ggplot(corr_q_vals_summ %>% filter(probability=="90-10"),   
       aes(x=trial_within_condition, y=m, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ cond) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.9, color="blue", size=1.5, linetype="longdash") +
  geom_hline(yintercept=-.1, color="red", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("Q-value") +
  ggtitle("Correct") +
  xlab("Stim iteration") + ylim(-1, 1)
b <- 
  ggplot(incorr_q_vals_summ %>% filter(probability=="90-10"),   
       aes(x=trial_within_condition, y=m, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ cond) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.1, color="blue", size=1.5, linetype="longdash") +
  geom_hline(yintercept=-.9, color="red", size=1.5, linetype="longdash") +
  ga + ap + tp + ft + tol +
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ggtitle("Incorrect") +
  ylab("") +
  xlab("Stim iteration") + ylim(-1, 1)
c <- 
   ggplot(corr_q_vals_summ %>% filter(probability=="90-10"),   
       aes(x=trial_within_condition, y=diff, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ cond) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.8, color="blue", size=1.5, linetype="longdash") +
  geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("") +
  ggtitle("Difference") +
  xlab("Stim iteration") + ylim(-1, 1)

# For legend  
# ggplot(corr_q_vals_summ %>% filter(probability=="90-10"),   
#        aes(x=trial_within_condition, y=diff, color=valence)) +
#   geom_hline(yintercept=0, color="black", size=2) + 
#   geom_line(size=3) + facet_wrap( ~ cond) + 
#   scale_color_manual(values=c("blue", "red")) +
#   geom_hline(yintercept=.8, color="blue", size=1.5, linetype="longdash") +
#   geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
#   ga + ap +   ft + 
#   theme(legend.text = element_text(size = 30),
#                legend.title = element_blank(),
#                legend.key.size = unit(2.5, 'lines')) +
#   theme(plot.title = element_text(size = 30, hjust = .5)) +
#   ylab("") +
#   ggtitle("Difference") +
#   xlab("Stim iteration") + ylim(-1, 1)
d <- ggplot(sim_perf_summ %>% filter(probability=="90-10"),   
       aes(x=trial_within_condition, y=m, color=valence)) +
  geom_hline(yintercept=0.5, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ cond) + 
  scale_color_manual(values=c("blue", "red")) +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("Proportion correct") +
  ggtitle("Learning") +
  xlab("Stim iteration") 
sim_test_perf_summ <- 
  s1_test_sim_m1_eb %>% filter(ID==25) %>% 
  group_by(cond, probability, valence) %>% summarize(m=mean(sit_sim_corrs))
## `summarise()` has grouped output by 'cond', 'probability'. You can override
## using the `.groups` argument.
sim_test_perf_summ$valence <- factor(sim_test_perf_summ$valence, levels=c("reward", "punishment"))

sim_test_perf_summ[sim_test_perf_summ$cond=="cog", "cond"] <- "cognitive"
#sim_test_perf_summ %>% filter(probability=="90-10")
e <- ggplot(sim_test_perf_summ %>% filter(probability=="90-10"),   
       aes(x=valence, y=m, color=valence)) + 
  geom_hline(yintercept=0.5, color="black", size=2) +
  geom_bar(stat="identity", fill="white", size=3) + 
  facet_wrap(~cond) +
  
  scale_color_manual(values=c("blue", "red")) +
  ga + ap + tol + ft +
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("Proportion correct") +
  ggtitle("Test") +
  
  xlab("") + theme(axis.ticks.x=element_blank(), axis.text.x=element_blank())#+ ylim(.48, 1) +
q_vals_comb <- a + b + c
q_vals_comb

sim_perf_comb <- d + e
sim_perf_comb

#ggsave("../paper/figs/fig_parts/q_vals_plot.png", q_vals_comb, height=5, width=16, dpi=700)
#ggsave("../paper/figs/fig_parts/sim_perf_comb.png", sim_perf_comb, height=6, width=14, dpi=700)

Supplementary Figure 7

Plotting Q-values by valence for m3 vs. m4 (m6 vs. m7) to understand how the introduction of pessimistic priors affects performance

# Pessimistic  
m3_sims_s1 <-  
  rm(sp,  "SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDecayToPessPrior82512.csv")
# Naive  
m4_sims_s1 <-  
  rm(sp, "SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDecayTo0Inits36394.csv")

M6: Q nothing varied in paper

pess_priors_qv_summs <- m3_sims_s1 %>% filter(probability=="90-10") %>% 
  group_by(valence, trial_within_condition) %>% 
  summarize(m_corr=mean(sim_correct_Q_vals), m_incorr=mean(sim_incorrect_Q_vals)) %>% 
  mutate(diff=m_corr-m_incorr)
## `summarise()` has grouped output by 'valence'. You can override using the
## `.groups` argument.
pess_priors_qv_summs$valence <- factor(pess_priors_qv_summs$valence, levels=c("reward", "punishment"))

M7: Q nothing varied, decay to 0

naive_priors_qv_summs <- m4_sims_s1 %>% filter(probability=="90-10") %>% 
  group_by(valence, trial_within_condition) %>% 
  summarize(m_corr=mean(sim_correct_Q_vals), m_incorr=mean(sim_incorrect_Q_vals)) %>% 
  mutate(diff=m_corr-m_incorr)
## `summarise()` has grouped output by 'valence'. You can override using the
## `.groups` argument.
naive_priors_qv_summs$valence <- factor(naive_priors_qv_summs$valence, levels=c("reward", "punishment"))

Break down into correct/incorrect plots

ggplot(pess_priors_qv_summs,   
       aes(x=trial_within_condition, y=m_corr, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ valence) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
  # geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("") +
  ggtitle("Correct") +
  xlab("Stim iteration") + ylim(-1, 1)

ggplot(naive_priors_qv_summs,   
       aes(x=trial_within_condition, y=m_corr, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ valence) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
  # geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("") +
  ggtitle("Correct") +
  xlab("Stim iteration") + ylim(-1, 1)

ggplot(pess_priors_qv_summs,   
       aes(x=trial_within_condition, y=m_incorr, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ valence) + 
  scale_color_manual(values=c("blue", "red")) +
  #geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
  # geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("") +
  ggtitle("Incorrect") +
  xlab("Stim iteration") + ylim(-1, 1)

ggplot(naive_priors_qv_summs,   
       aes(x=trial_within_condition, y=m_incorr, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ valence) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
  # geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("") +
  ggtitle("Incorrect") +
  xlab("Stim iteration") + ylim(-1, 1)

pess_plot <- ggplot(pess_priors_qv_summs,   
       aes(x=trial_within_condition, y=diff, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ valence) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("Q-value Difference: \n Correct vs. Incorrect") +
  ggtitle("Pessimistic Priors") +
  xlab("Stim iteration") + ylim(0, 1)
naive_plot <- ggplot(naive_priors_qv_summs,   
       aes(x=trial_within_condition, y=diff, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ valence) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("") +
  ggtitle("Naive Priors (All Zero)") +
  xlab("Stim iteration") + ylim(0, 1)
simple_qv_plots_comb <- pess_plot + naive_plot
simple_qv_plots_comb 

No powerpoint — loaded directly into supplemental

#ggsave("../paper/figs/fig_parts/simple_qv_plots_comb.png", simple_qv_plots_comb, height=7, width=14, dpi=700)

PST/Stimulus valuation phase analyses

Used rew history because Q-values are for specific Q(s,a)s whereas selection in the PST phase is between stimuli ie. Q(s, :)

# Take just the trials within condition   
within_cond_pst_s1 <- s1_pst %>% filter(test_condition == "cogn_cogn" | test_condition == "overt_overt")
within_cond_pst_s1$test_condition <- factor(within_cond_pst_s1$test_condition, levels=c("overt_overt", "cogn_cogn"))
contrasts(within_cond_pst_s1$test_condition) <- c(-.5, .5)
head(within_cond_pst_s1$test_condition)
## [1] overt_overt overt_overt overt_overt overt_overt overt_overt overt_overt
## attr(,"contrasts")
##             [,1]
## overt_overt -0.5
## cogn_cogn    0.5
## Levels: overt_overt cogn_cogn
# Both singular  
# summary(m0_s1 <- glmer(resp_num ~  scale(left_min_right) + ( scale(left_min_right)|ID),
#           data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))
# 
# summary(m0_s1 <- glmer(resp_num ~  scale(left_min_right)*test_condition + (scale(left_min_right) |ID),
#           data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))

Using model without random intercepts because not singular

summary(m0_s1 <- glmer(resp_num ~  scale(left_min_right) + (0 + scale(left_min_right)|ID),
          data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: resp_num ~ scale(left_min_right) + (0 + scale(left_min_right) |  
##     ID)
##    Data: within_cond_pst_s1
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   6456.5   6476.6  -3225.2   6450.5     5983 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3395 -0.7468  0.1149  0.7468  4.5130 
## 
## Random effects:
##  Groups Name                  Variance Std.Dev.
##  ID     scale(left_min_right) 1.538    1.24    
## Number of obs: 5986, groups:  ID, 125
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.11210    0.03073   3.648 0.000264 ***
## scale(left_min_right)  1.41202    0.12083  11.686  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(lft_m_) 0.015
summary(m1_s1 <- glmer(resp_num ~  scale(left_min_right)*test_condition + (0 + scale(left_min_right) |ID),
          data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## resp_num ~ scale(left_min_right) * test_condition + (0 + scale(left_min_right) |  
##     ID)
##    Data: within_cond_pst_s1
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   6431.2   6464.7  -3210.6   6421.2     5981 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.1380 -0.7457  0.1136  0.7448  4.1139 
## 
## Random effects:
##  Groups Name                  Variance Std.Dev.
##  ID     scale(left_min_right) 1.555    1.247   
## Number of obs: 5986, groups:  ID, 125
## 
## Fixed effects:
##                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                            0.11079    0.03085   3.592 0.000328 ***
## scale(left_min_right)                  1.42529    0.12149  11.731  < 2e-16 ***
## test_condition1                        0.09218    0.06167   1.495 0.134952    
## scale(left_min_right):test_condition1 -0.36515    0.07029  -5.195 2.05e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sc(__) tst_c1
## scl(lft_m_)  0.015              
## test_cndtn1 -0.041  0.005       
## scl(l__):_1  0.003 -0.035  0.025
# Take just the trials within condition   
within_cond_pst_s2 <- s2_pst %>% filter(test_condition == "cognitive_cognitive" | test_condition == "overt_overt")
within_cond_pst_s2$test_condition <- 
  factor(within_cond_pst_s2$test_condition, levels=c("overt_overt", "cognitive_cognitive"))
contrasts(within_cond_pst_s2$test_condition) <- c(-.5, .5)
head(within_cond_pst_s2$test_condition)
## [1] cognitive_cognitive cognitive_cognitive cognitive_cognitive
## [4] cognitive_cognitive cognitive_cognitive cognitive_cognitive
## attr(,"contrasts")
##                     [,1]
## overt_overt         -0.5
## cognitive_cognitive  0.5
## Levels: overt_overt cognitive_cognitive

Use the comparable model

summary(m0_s2 <- glmer(resp_num ~  scale(left_min_right) + (0 + scale(left_min_right)|ID),
          data=within_cond_pst_s2, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: resp_num ~ scale(left_min_right) + (0 + scale(left_min_right) |  
##     ID)
##    Data: within_cond_pst_s2
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   7161.3   7181.7  -3577.7   7155.3     6616 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.1402 -0.7899  0.1255  0.7803  4.7614 
## 
## Random effects:
##  Groups Name                  Variance Std.Dev.
##  ID     scale(left_min_right) 1.507    1.228   
## Number of obs: 6619, groups:  ID, 138
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.11126    0.02914   3.818 0.000134 ***
## scale(left_min_right)  1.35153    0.11346  11.912  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(lft_m_) 0.014
summary(m1_s2 <- glmer(resp_num ~  scale(left_min_right)*test_condition + (0 + scale(left_min_right) |ID),
          data=within_cond_pst_s2, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## resp_num ~ scale(left_min_right) * test_condition + (0 + scale(left_min_right) |  
##     ID)
##    Data: within_cond_pst_s2
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   7163.5   7197.5  -3576.8   7153.5     6614 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.0975 -0.7881  0.1258  0.7749  4.9593 
## 
## Random effects:
##  Groups Name                  Variance Std.Dev.
##  ID     scale(left_min_right) 1.509    1.228   
## Number of obs: 6619, groups:  ID, 138
## 
## Fixed effects:
##                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                            0.11101    0.02914   3.809  0.00014 ***
## scale(left_min_right)                  1.35180    0.11350  11.910  < 2e-16 ***
## test_condition1                        0.06205    0.05826   1.065  0.28689    
## scale(left_min_right):test_condition1 -0.05321    0.06672  -0.797  0.42519    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sc(__) tst_c1
## scl(lft_m_)  0.014              
## test_cndtn1 -0.008  0.004       
## scl(l__):_1  0.006 -0.002  0.027
# Singular 
# summary(m1_s2 <- glmer(resp_num ~  scale(left_min_right)*test_condition + (0 + scale(left_min_right)*test_condition |ID),
#           data=within_cond_pst_s2, family="binomial", control = glmerControl(optimizer = "bobyqa")))
p_s1 <- 
  sjPlot::plot_model(m1_s1, type="pred", terms = c("left_min_right [all]", "test_condition"))


p_fin_1 <- p_s1 + ap + tp + ga + 
  xlab("Left minus right reward history") + ylab("") + ggtitle("") + 
  theme(plot.title = element_text(hjust=0)) +
  ylab("Chance of picking left") +
   xlab("Left minus right reward history") +
  ggtitle("Predicted probabilities") + 
  scale_color_manual(values=c("orange", "purple"), labels=c("overt-overt", "cognitive-cognitive")) + tol
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
p_s2 <- 
  sjPlot::plot_model(m1_s2, type="pred", terms = c("left_min_right [all]", "test_condition"))


p_fin_s2 <- p_s2 + ap + tp + ga + 
  xlab("Left minus right reward history") + ylab("") + ggtitle("") + 
  theme(plot.title = element_text(hjust=0)) +
  ylab("Chance of picking left") +
   xlab("Left minus right reward history") +
  #ggtitle("Predicted probabilities") + 
  scale_color_manual(values=c("orange", "purple"), labels=c("overt-overt", "cognitive-cognitive")) +
  theme(legend.text = element_text(size = 18),
               legend.title = element_blank(),
               legend.key.size = unit(1.3, 'lines'))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
p_fin_1 + p_fin_s2

ranef(m1_s2)
## $ID
##     scale(left_min_right)
## 1           -7.638427e-02
## 2           -1.911542e-01
## 3           -1.535659e+00
## 4            1.398870e-01
## 5            1.763300e+00
## 6           -1.465558e+00
## 7           -5.962328e-01
## 8            8.459447e-01
## 9           -2.737737e-01
## 10           8.023781e-01
## 11          -1.789991e+00
## 12          -3.571217e-01
## 13          -1.554381e+00
## 14           6.793301e-01
## 15           3.156066e-01
## 16          -1.700867e-01
## 17          -1.283673e+00
## 18          -9.527213e-01
## 19          -7.157713e-01
## 20          -4.166149e-01
## 21          -1.332563e+00
## 22          -9.936151e-01
## 23           4.158967e-01
## 24          -1.262534e+00
## 25          -8.249032e-01
## 26           1.347361e+00
## 27          -1.197001e+00
## 28           1.846921e+00
## 29          -1.252705e+00
## 30           6.828607e-01
## 31          -1.364753e+00
## 32          -1.792755e+00
## 33           8.582025e-01
## 34          -1.775080e-01
## 35          -5.372453e-02
## 36           6.800819e-01
## 37           1.997015e-01
## 38          -5.007704e-01
## 39           1.040108e+00
## 40           2.002604e+00
## 41          -1.831436e+00
## 42           6.286168e-01
## 43           1.583285e+00
## 44          -7.409308e-01
## 45           1.068599e+00
## 46          -3.780431e-01
## 47           6.087241e-01
## 48           2.624062e-01
## 49           1.076875e+00
## 50          -1.681045e+00
## 51           7.740095e-02
## 52           1.276172e-01
## 53          -8.255441e-01
## 54           1.501652e+00
## 55          -1.198909e+00
## 56           1.231464e+00
## 57          -7.374282e-01
## 58          -1.044010e+00
## 59           1.330937e+00
## 60           9.092997e-01
## 61           6.311263e-03
## 62           1.141079e+00
## 63          -6.813984e-01
## 64          -6.699732e-01
## 65          -1.327898e+00
## 66          -1.787861e-01
## 67          -1.396548e-01
## 68          -3.535874e-01
## 69           1.867799e-01
## 70          -1.221560e+00
## 71          -5.083901e-01
## 72           9.698690e-01
## 73           6.657423e-01
## 74          -7.371641e-01
## 75          -1.415828e+00
## 76          -1.264396e+00
## 77           1.724080e+00
## 78          -1.221658e+00
## 79           1.657571e+00
## 80          -1.628377e+00
## 81           2.921264e-01
## 82          -2.411008e-01
## 83          -9.298301e-05
## 84           1.107018e+00
## 85          -9.797200e-02
## 86          -2.256282e+00
## 87           1.880165e+00
## 88          -1.238028e+00
## 89          -2.559056e+00
## 90          -1.093334e+00
## 91          -7.953529e-01
## 92          -4.959134e-01
## 93           1.412486e+00
## 94          -9.225457e-01
## 95           1.014909e+00
## 96          -8.901824e-01
## 97           2.365885e-01
## 98          -1.210926e+00
## 99           3.047484e-02
## 100          3.144635e-01
## 101         -1.402683e+00
## 102         -8.816298e-01
## 103          1.511560e+00
## 104         -9.874832e-01
## 105          2.173403e-01
## 106         -1.052314e+00
## 107          1.724758e+00
## 108          1.839174e+00
## 109         -1.697650e+00
## 110         -9.228177e-01
## 111          1.453432e+00
## 112          8.494656e-01
## 113          2.214903e+00
## 114         -4.507024e-01
## 115         -3.275397e-01
## 116         -1.026534e+00
## 117          7.187065e-01
## 118          2.520358e-01
## 119          2.038057e+00
## 120         -1.233438e+00
## 121          1.897789e+00
## 122         -1.107078e+00
## 123          1.644766e+00
## 124          9.238578e-01
## 125          9.303096e-02
## 126         -2.417370e-01
## 127          1.273118e+00
## 128          6.328269e-01
## 129          1.196723e+00
## 130         -1.583174e-01
## 131         -1.289770e+00
## 132         -5.052138e-01
## 133          7.200495e-01
## 134          9.689824e-01
## 135          1.159983e+00
## 136         -1.270916e+00
## 137         -2.006149e-01
## 138          1.017842e+00
## 
## with conditional variances for "ID"

Mixed test phase

unique(s1_pst$test_condition)
## [1] "overt_cogn"  "overt_overt" "cogn_cogn"
mix_s1 <- s1_pst %>% filter(test_condition == "overt_cogn")
# If cognitive is on the left and they chose left then they chose cognitive  
mix_s1[mix_s1$left_training_cond=="cognTraining", "chose_cognitive"] <- 
  if_else(mix_s1[mix_s1$left_training_cond=="cognTraining", "response"]=="left", 1, 0)

# If cognitive is on the right and they chose right then they chose cognitive
mix_s1[mix_s1$right_training_cond=="cognTraining", "chose_cognitive"] <- 
  if_else(mix_s1[mix_s1$right_training_cond=="cognTraining", "response"]=="right", 1, 0)
table(mix_s1$chose_cognitive)[2]/sum(table(mix_s1$chose_cognitive))
##         1 
## 0.4751506
unique(s2_pst$test_condition)
## [1] "overt_cognitive"     "cognitive_cognitive" "overt_overt"
mix_s2 <- s2_pst %>% filter(test_condition == "overt_cognitive")
# If cognitive is on the left and they chose left then they chose cognitive  
mix_s2[mix_s2$left_training_cond=="cognitive", "chose_cognitive"] <- 
  if_else(mix_s2[mix_s2$left_training_cond=="cognitive", "response"]=="left", 1, 0)

# If cognitive is on the right and they chose right then they chose cognitive
mix_s2[mix_s2$right_training_cond=="cognitive", "chose_cognitive"] <- 
  if_else(mix_s2[mix_s2$right_training_cond=="cognitive", "response"]=="right", 1, 0)
table(mix_s2$chose_cognitive)[2]/sum(table(mix_s2$chose_cognitive))
##         1 
## 0.4231989
mixed_summs_s1 <- mix_s1 %>% group_by(choice_type_plotting) %>% 
  summarize(m=mean(chose_cognitive), n())
mixed_summs_s1
## # A tibble: 4 × 3
##   choice_type_plotting     m `n()`
##   <chr>                <dbl> <int>
## 1 "freq-P\nfreq-P"     0.490   997
## 2 "freq-R\nfreq-R"     0.428   997
## 3 "infreq-P\ninfreq-P" 0.553   992
## 4 "infreq-R\ninfreq-R" 0.429   998
mixed_summs_s2 <- mix_s2 %>% group_by(choice_type_plotting) %>% 
  summarize(m=mean(chose_cognitive), n())
mixed_summs_s2
## # A tibble: 4 × 3
##   choice_type_plotting     m `n()`
##   <chr>                <dbl> <int>
## 1 "freq-P\nfreq-P"     0.408  1104
## 2 "freq-R\nfreq-R"     0.393  1104
## 3 "infreq-P\ninfreq-P" 0.447  1102
## 4 "infreq-R\ninfreq-R" 0.445  1104
summary(choice_effect_s1 <- glmer(chose_cognitive ~  1 + ( 1 |ID),
          data=mix_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: chose_cognitive ~ 1 + (1 | ID)
##    Data: mix_s1
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   5269.0   5281.6  -2632.5   5265.0     3982 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4732 -0.8937 -0.4973  0.9600  2.3488 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 0.5492   0.7411  
## Number of obs: 3984, groups:  ID, 125
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.10926    0.07443  -1.468    0.142
summary(choice_effect_s2 <- glmer(chose_cognitive ~  1 + ( 1|ID),
          data=mix_s2, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: chose_cognitive ~ 1 + (1 | ID)
##    Data: mix_s2
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   5829.4   5842.2  -2912.7   5825.4     4412 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5760 -0.7917 -0.5754  1.0417  2.2498 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 0.3916   0.6258  
## Number of obs: 4414, groups:  ID, 138
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.33769    0.06215  -5.433 5.53e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3/19/24 - ethnicity as per Ethnicity.identified var from Prolific from shared drive (shared drive doc has other demographic tallies)

# Study 1
s1n <- 125
7/s1n # Asian
## [1] 0.056
8/s1n # Black
## [1] 0.064
5/s1n # Mixed
## [1] 0.04
3/s1n # Other
## [1] 0.024
100/s1n # White
## [1] 0.8
2/s1n # Data expired
## [1] 0.016
# Study 2
s2n <- 4+5+6+13+5+105
s2n==138
## [1] TRUE
4/s2n # Asian
## [1] 0.02898551
5/s2n # Black
## [1] 0.03623188
13/s2n # Mixed
## [1] 0.0942029
5/s2n # Other
## [1] 0.03623188
105/s2n # White
## [1] 0.7608696
6/s2n # Data expired
## [1] 0.04347826

Questionnaires

s1_qs <- read.csv("../data/cleaned_files/s1_qs_deidentified.csv")
s2_qs <- read.csv("../data/cleaned_files/s2_qs_deidentified.csv")
s1_demogs <- read.csv("../data/cleaned_files/s1_demogs_deIDed.csv")
s2_demogs <- read.csv("../data/cleaned_files/s2_demogs_deIDed.csv")
# names(s1_qs)
# names(s2_qs)

Exclude pts who answered both catch questions incorrectly

S1 all correct so no exclusions

table(s1_qs$Q6_20) # All correct
## 
## A little bit 
##          123
table(s1_qs$Q7_19) # All correct
## 
## Very false for me 
##               123

4/4 — now removing the pt who skipped all items here (leaving n=122 completers). This mayslightly influenced PTQ and RRS imputations; statistic now corrected to p’s > .07 (originally .10)

s1_qs <- s1_qs[-66, ]
table(s2_qs$Q6_20) # All correct
## 
## A little bit 
##          137
table(s2_qs$Q7_19) # One pt incorrect but they answered the other catch item correctly, so retained (given recs in lit that excluding on a single catch item can bias results)
## 
## Very false for me  Very true for me 
##               136                 1
#s2_qs %>% filter(Q7_19 == "Very true for me")

From notes:
Q2 is the PTQ
Q6 is the MASQ — 23 items, catch is 20 (confirmed on Qualtrics) - but permission issue explained in revision hence not reported
Q7 is the BIS/BAS — 25 items, catch is 19 (confirmed on Qualtrics)
- some are not part of the scale
Q8 is the RRS — 10 items, no catch

S1 questionnaires

Get and recode PTQ and RRS

s1_ptq <- s1_qs %>% select(contains("Q2_")) %>% setNames(paste0("PTQ_", 1:15))
#s1_ptq$ID <- s1_qs$ID 4/4 removing to make sure not in imputation 
s1_ptq_copy <- s1_ptq
s1_ptq_copy$ID <- s1_qs$ID
s1_ptq[s1_ptq=="Never"] <- 0
s1_ptq[s1_ptq=="Rarely"] <- 1
s1_ptq[s1_ptq=="Sometimes"] <- 2 
s1_ptq[s1_ptq=="Often"] <- 3
s1_ptq[s1_ptq=="Almost always"] <- 4
s1_rrs <- 
  s1_qs %>% select(contains("Q8_")) %>% setNames(paste0("RRS_", 1:10))
#s1_rrs$ID <- s1_qs$ID 4/4 removing to make sure not in imputation 
s1_rrs_copy <- s1_rrs
s1_rrs_copy$ID <- s1_qs$ID
s1_rrs[s1_rrs=="Almost never"] <- 1
s1_rrs[s1_rrs=="Sometimes"] <- 2
s1_rrs[s1_rrs=="Often"] <- 3
s1_rrs[s1_rrs=="Almost always"] <- 4

For revision: BIS/BAS

Exclude catch but keep filler items for now

#s1_qs %>% select(contains("Q7_")) %>% select("Q7_19") # confirmed this is catch 
s1_bis_bas <- s1_qs %>% select(contains("Q7_")) %>% select(-"Q7_19") 

“Items other than 2 and 22 are reverse coded”

The two BIS items that are forward coded (1= very true, 4=very false) — Q7_23 corresponds to item 22 in actual scale because of the catch item

fwd_code_bis_bas <- s1_bis_bas %>% select("Q7_2", "Q7_23")
fwd_code_bis_bas[fwd_code_bis_bas=="Very true for me"] <- 1
fwd_code_bis_bas[fwd_code_bis_bas=="Somewhat true for me"] <- 2
fwd_code_bis_bas[fwd_code_bis_bas=="Somewhat false for me"] <- 3
fwd_code_bis_bas[fwd_code_bis_bas=="Very false for me"] <- 4

Sanity check that positively correlate

cor.test(as.numeric(fwd_code_bis_bas$Q7_2), as.numeric(fwd_code_bis_bas$Q7_23))
## 
##  Pearson's product-moment correlation
## 
## data:  as.numeric(fwd_code_bis_bas$Q7_2) and as.numeric(fwd_code_bis_bas$Q7_23)
## t = 8.8719, df = 120, p-value = 8.218e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5085008 0.7259186
## sample estimates:
##       cor 
## 0.6293718

The rest are reverse coded

rev_code_bis_bas <- s1_bis_bas %>% select(-c("Q7_2", "Q7_23"))
rev_code_bis_bas[rev_code_bis_bas=="Very true for me"] <- 4
rev_code_bis_bas[rev_code_bis_bas=="Somewhat true for me"] <- 3
rev_code_bis_bas[rev_code_bis_bas=="Somewhat false for me"] <- 2
rev_code_bis_bas[rev_code_bis_bas=="Very false for me"] <- 1
s1_bis_bas <- data.frame(fwd_code_bis_bas, rev_code_bis_bas)
s1_bis_bas_num <- foreach (i = 1:ncol(s1_bis_bas)) %do% {
  data.frame(as.numeric(unlist(s1_bis_bas[i]))) %>% setNames(names(s1_bis_bas[i]))
}%>% bind_cols()

Sanity check BIS fwd and reverse code items postiively correlate

cor.test(s1_bis_bas_num$Q7_2, s1_bis_bas_num$Q7_23)
## 
##  Pearson's product-moment correlation
## 
## data:  s1_bis_bas_num$Q7_2 and s1_bis_bas_num$Q7_23
## t = 8.8719, df = 120, p-value = 8.218e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5085008 0.7259186
## sample estimates:
##       cor 
## 0.6293718
cor.test(s1_bis_bas_num$Q7_2, s1_bis_bas_num$Q7_8)
## 
##  Pearson's product-moment correlation
## 
## data:  s1_bis_bas_num$Q7_2 and s1_bis_bas_num$Q7_8
## t = 6.3919, df = 120, p-value = 3.264e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3583137 0.6256823
## sample estimates:
##       cor 
## 0.5039744
cor.test(s1_bis_bas_num$Q7_23, s1_bis_bas_num$Q7_8)
## 
##  Pearson's product-moment correlation
## 
## data:  s1_bis_bas_num$Q7_23 and s1_bis_bas_num$Q7_8
## t = 4.8078, df = 120, p-value = 4.47e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2413693 0.5409989
## sample estimates:
##       cor 
## 0.4018868
# Item 24 because shifted by 1 due to catch
cor.test(s1_bis_bas_num$Q7_2, s1_bis_bas_num$Q7_25)
## 
##  Pearson's product-moment correlation
## 
## data:  s1_bis_bas_num$Q7_2 and s1_bis_bas_num$Q7_25
## t = 8.6941, df = 120, p-value = 2.152e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4990482 0.7198710
## sample estimates:
##       cor 
## 0.6216608
cor.test(s1_bis_bas_num$Q7_23, s1_bis_bas_num$Q7_25)
## 
##  Pearson's product-moment correlation
## 
## data:  s1_bis_bas_num$Q7_23 and s1_bis_bas_num$Q7_25
## t = 5.8332, df = 120, p-value = 4.707e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3188964 0.5978260
## sample estimates:
##       cor 
## 0.4700135
cor.test(s1_bis_bas_num$Q7_2, s1_bis_bas_num$Q7_16)
## 
##  Pearson's product-moment correlation
## 
## data:  s1_bis_bas_num$Q7_2 and s1_bis_bas_num$Q7_16
## t = 6.2646, df = 120, p-value = 6.056e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3495145 0.6195225
## sample estimates:
##       cor 
## 0.4964321
cor.test(s1_bis_bas_num$Q7_23, s1_bis_bas_num$Q7_16)
## 
##  Pearson's product-moment correlation
## 
## data:  s1_bis_bas_num$Q7_23 and s1_bis_bas_num$Q7_16
## t = 6.0606, df = 120, p-value = 1.611e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3351901 0.6094231
## sample estimates:
##       cor 
## 0.4841061

Now remove filler items — all before catch so all items correspond to actual num in scale

s1_bis_bas <- s1_bis_bas_num %>% select(-c(
  "Q7_1",
  "Q7_6",
  "Q7_11",
  "Q7_17") 
)

BIS

s1_bis <- s1_bis_bas %>% select(
  # Before catch item 
  c("Q7_2", "Q7_8", "Q7_13", "Q7_16"),
  # After catch so one shifted
  c("Q7_20", "Q7_23", "Q7_25"))

Just 1 missing so impute

length(which(is.na(s1_bis[, 1:7])))
## [1] 1
length(which(is.na(s1_bis[, 1:7])))/(dim(s1_bis[, 1:7])[1]*dim(s1_bis[, 1:7])[2]) 
## [1] 0.00117096
s1_bis_final <- foreach (i = 1:nrow(s1_bis)) %do% {
  if (any(is.na(s1_bis[i, ]))) {
    s1_bis[i, is.na(s1_bis[i, ])] <- mean(unlist(s1_bis[i, ]), na.rm=TRUE)
  }
s1_bis[i, ]
}%>% bind_rows()

Sanity check all positively correlated

cor(s1_bis_final)
##            Q7_2      Q7_8     Q7_13     Q7_16     Q7_20     Q7_23     Q7_25
## Q7_2  1.0000000 0.5039744 0.5156167 0.4964321 0.5391917 0.6293718 0.6216608
## Q7_8  0.5039744 1.0000000 0.6775552 0.5204724 0.5118818 0.4018868 0.5865127
## Q7_13 0.5156167 0.6775552 1.0000000 0.5338191 0.5322763 0.4445993 0.4822497
## Q7_16 0.4964321 0.5204724 0.5338191 1.0000000 0.5482658 0.4841061 0.5369988
## Q7_20 0.5391917 0.5118818 0.5322763 0.5482658 1.0000000 0.3565339 0.6321968
## Q7_23 0.6293718 0.4018868 0.4445993 0.4841061 0.3565339 1.0000000 0.4700135
## Q7_25 0.6216608 0.5865127 0.4822497 0.5369988 0.6321968 0.4700135 1.0000000

And high alpha

psych::alpha(s1_bis_final)
## 
## Reliability analysis   
## Call: psych::alpha(x = s1_bis_final)
## 
##   raw_alpha std.alpha G6(smc) average_r S/N   ase mean  sd median_r
##       0.88      0.89    0.88      0.53 7.7 0.016    3 0.7     0.52
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt     0.85  0.88  0.91
## Duhachek  0.85  0.88  0.91
## 
##  Reliability if an item is dropped:
##       raw_alpha std.alpha G6(smc) average_r S/N alpha se  var.r med.r
## Q7_2       0.86      0.86    0.86      0.51 6.4    0.020 0.0067  0.52
## Q7_8       0.86      0.87    0.86      0.52 6.5    0.019 0.0053  0.53
## Q7_13      0.86      0.87    0.86      0.52 6.6    0.019 0.0061  0.52
## Q7_16      0.87      0.87    0.87      0.53 6.7    0.019 0.0082  0.52
## Q7_20      0.87      0.87    0.86      0.53 6.7    0.019 0.0055  0.52
## Q7_23      0.88      0.88    0.87      0.55 7.3    0.017 0.0031  0.53
## Q7_25      0.86      0.86    0.86      0.51 6.3    0.019 0.0062  0.52
## 
##  Item statistics 
##         n raw.r std.r r.cor r.drop mean   sd
## Q7_2  122  0.80  0.80  0.76   0.71  3.0 0.90
## Q7_8  122  0.79  0.78  0.74   0.69  2.8 1.01
## Q7_13 122  0.79  0.78  0.74   0.69  2.9 0.99
## Q7_16 122  0.76  0.76  0.71   0.67  2.8 0.88
## Q7_20 122  0.75  0.76  0.72   0.66  3.3 0.81
## Q7_23 122  0.71  0.70  0.64   0.59  2.8 1.00
## Q7_25 122  0.79  0.80  0.77   0.72  3.2 0.83
## 
## Non missing response frequency for each item
##          1    2 2.66666666666667    3    4 miss
## Q7_2  0.07 0.18             0.00 0.39 0.36    0
## Q7_8  0.13 0.25             0.00 0.33 0.29    0
## Q7_13 0.11 0.20             0.00 0.34 0.34    0
## Q7_16 0.07 0.30             0.00 0.41 0.22    0
## Q7_20 0.06 0.06             0.01 0.45 0.43    0
## Q7_23 0.12 0.25             0.00 0.34 0.29    0
## Q7_25 0.04 0.14             0.00 0.41 0.41    0
s1_bis_final$bis_sum <- rowSums(s1_bis_final)
s1_bis_final$ID <- s1_qs$ID

The remaining items are BAS

s1_bas <- s1_bis_bas %>% select(
  
  -c("Q7_2", "Q7_8", "Q7_13", "Q7_16",
  
  "Q7_20", "Q7_23", "Q7_25"))

Just 3 missing so impute

length(which(is.na(s1_bas[, 1:ncol(s1_bas)])))
## [1] 3
length(which(is.na(s1_bas[, 1:ncol(s1_bas)])))/(dim(s1_bas[, 1:ncol(s1_bas)])[1]*dim(s1_bas[, 1:ncol(s1_bas)])[2]) 
## [1] 0.001891551
s1_bas_final <- foreach (i = 1:nrow(s1_bas)) %do% {
  if (any(is.na(s1_bas[i, ]))) {
    s1_bas[i, is.na(s1_bas[i, ])] <- mean(unlist(s1_bas[i, ]), na.rm=TRUE)
  }
s1_bas[i, ]
}%>% bind_rows()

Sanity check all positively correlated

cor(s1_bas_final)
##            Q7_3      Q7_4      Q7_5      Q7_7      Q7_9     Q7_10     Q7_12
## Q7_3  1.0000000 0.3505651 0.5194718 0.3841874 0.6711250 0.5025063 0.6169502
## Q7_4  0.3505651 1.0000000 0.2879864 0.3390419 0.2583916 0.1969002 0.3028181
## Q7_5  0.5194718 0.2879864 1.0000000 0.4144983 0.4738954 0.6536637 0.5278935
## Q7_7  0.3841874 0.3390419 0.4144983 1.0000000 0.3083571 0.3054471 0.2941940
## Q7_9  0.6711250 0.2583916 0.4738954 0.3083571 1.0000000 0.4900888 0.7280516
## Q7_10 0.5025063 0.1969002 0.6536637 0.3054471 0.4900888 1.0000000 0.3838347
## Q7_12 0.6169502 0.3028181 0.5278935 0.2941940 0.7280516 0.3838347 1.0000000
## Q7_14 0.5878563 0.3752357 0.5274854 0.6148655 0.5484696 0.3908887 0.6012681
## Q7_15 0.2808272 0.1617566 0.4120492 0.2799688 0.3525170 0.5021741 0.3237509
## Q7_18 0.2918551 0.1745012 0.1677867 0.5125142 0.2076797 0.1884286 0.2343946
## Q7_21 0.5414331 0.1512069 0.5872976 0.3294216 0.4492972 0.5246798 0.3926814
## Q7_22 0.5064174 0.2087005 0.3996587 0.1100529 0.6215593 0.3941751 0.5789322
## Q7_24 0.3864836 0.2818568 0.3438651 0.5351943 0.2465222 0.1690507 0.3441865
##           Q7_14     Q7_15     Q7_18     Q7_21     Q7_22     Q7_24
## Q7_3  0.5878563 0.2808272 0.2918551 0.5414331 0.5064174 0.3864836
## Q7_4  0.3752357 0.1617566 0.1745012 0.1512069 0.2087005 0.2818568
## Q7_5  0.5274854 0.4120492 0.1677867 0.5872976 0.3996587 0.3438651
## Q7_7  0.6148655 0.2799688 0.5125142 0.3294216 0.1100529 0.5351943
## Q7_9  0.5484696 0.3525170 0.2076797 0.4492972 0.6215593 0.2465222
## Q7_10 0.3908887 0.5021741 0.1884286 0.5246798 0.3941751 0.1690507
## Q7_12 0.6012681 0.3237509 0.2343946 0.3926814 0.5789322 0.3441865
## Q7_14 1.0000000 0.3246495 0.4751540 0.4521666 0.4152103 0.5021592
## Q7_15 0.3246495 1.0000000 0.2702716 0.5064189 0.3647177 0.1475499
## Q7_18 0.4751540 0.2702716 1.0000000 0.1857702 0.2346360 0.4633299
## Q7_21 0.4521666 0.5064189 0.1857702 1.0000000 0.4651129 0.3173153
## Q7_22 0.4152103 0.3647177 0.2346360 0.4651129 1.0000000 0.1378870
## Q7_24 0.5021592 0.1475499 0.4633299 0.3173153 0.1378870 1.0000000

And high alpha

psych::alpha(s1_bas_final) 
## 
## Reliability analysis   
## Call: psych::alpha(x = s1_bas_final)
## 
##   raw_alpha std.alpha G6(smc) average_r S/N   ase mean   sd median_r
##       0.89      0.89    0.92      0.39 8.2 0.014  2.9 0.52     0.38
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt     0.86  0.89  0.92
## Duhachek  0.86  0.89  0.92
## 
##  Reliability if an item is dropped:
##       raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## Q7_3       0.88      0.88    0.90      0.37 7.1    0.016 0.022  0.36
## Q7_4       0.89      0.89    0.92      0.41 8.3    0.015 0.021  0.41
## Q7_5       0.88      0.88    0.90      0.38 7.2    0.016 0.022  0.36
## Q7_7       0.89      0.88    0.91      0.39 7.6    0.015 0.023  0.39
## Q7_9       0.88      0.88    0.90      0.38 7.2    0.016 0.020  0.38
## Q7_10      0.88      0.88    0.91      0.38 7.5    0.015 0.022  0.37
## Q7_12      0.88      0.88    0.90      0.38 7.2    0.016 0.021  0.38
## Q7_14      0.88      0.87    0.90      0.37 7.0    0.016 0.022  0.35
## Q7_15      0.89      0.89    0.91      0.40 7.9    0.015 0.023  0.39
## Q7_18      0.89      0.89    0.91      0.40 8.2    0.014 0.021  0.39
## Q7_21      0.88      0.88    0.91      0.38 7.4    0.016 0.023  0.37
## Q7_22      0.88      0.88    0.91      0.39 7.6    0.015 0.021  0.38
## Q7_24      0.89      0.89    0.91      0.40 7.9    0.015 0.022  0.39
## 
##  Item statistics 
##         n raw.r std.r r.cor r.drop mean   sd
## Q7_3  122  0.78  0.78  0.76   0.72  2.9 0.81
## Q7_4  122  0.44  0.48  0.40   0.37  3.5 0.53
## Q7_5  122  0.75  0.74  0.72   0.68  3.0 0.83
## Q7_7  122  0.62  0.63  0.61   0.54  3.3 0.74
## Q7_9  122  0.75  0.74  0.74   0.69  2.8 0.79
## Q7_10 122  0.69  0.67  0.64   0.60  2.7 0.89
## Q7_12 122  0.74  0.74  0.73   0.68  2.8 0.76
## Q7_14 122  0.78  0.80  0.79   0.74  3.1 0.71
## Q7_15 122  0.59  0.58  0.53   0.50  2.2 0.82
## Q7_18 122  0.51  0.51  0.47   0.41  3.1 0.81
## Q7_21 122  0.71  0.69  0.67   0.63  2.6 0.85
## Q7_22 122  0.65  0.64  0.61   0.57  2.3 0.86
## Q7_24 122  0.55  0.57  0.53   0.47  3.4 0.70
## 
## Non missing response frequency for each item
##          1    2 2.91666666666667    3 3.41666666666667 3.58333333333333    4
## Q7_3  0.07 0.20             0.00 0.52             0.00             0.00 0.20
## Q7_4  0.00 0.02             0.00 0.47             0.01             0.00 0.51
## Q7_5  0.04 0.20             0.00 0.43             0.00             0.00 0.32
## Q7_7  0.03 0.07             0.00 0.47             0.00             0.00 0.43
## Q7_9  0.04 0.33             0.00 0.45             0.00             0.01 0.17
## Q7_10 0.10 0.31             0.00 0.41             0.00             0.00 0.18
## Q7_12 0.03 0.29             0.00 0.49             0.00             0.00 0.19
## Q7_14 0.02 0.15             0.00 0.54             0.00             0.00 0.30
## Q7_15 0.19 0.45             0.00 0.30             0.00             0.00 0.06
## Q7_18 0.05 0.14             0.01 0.48             0.00             0.00 0.32
## Q7_21 0.09 0.33             0.00 0.43             0.00             0.00 0.16
## Q7_22 0.18 0.41             0.00 0.33             0.00             0.00 0.08
## Q7_24 0.02 0.07             0.00 0.39             0.00             0.00 0.52
##       miss
## Q7_3     0
## Q7_4     0
## Q7_5     0
## Q7_7     0
## Q7_9     0
## Q7_10    0
## Q7_12    0
## Q7_14    0
## Q7_15    0
## Q7_18    0
## Q7_21    0
## Q7_22    0
## Q7_24    0
s1_bas_final$bas_sum <- rowSums(s1_bas_final)
hist(s1_bas_final$bas_sum, breaks=25)

s1_bas_final$ID <- s1_qs$ID

Modestly inversely related

ComparePars(s1_bas_final$bas_sum, s1_bis_final$bis_sum, use_identity_line = 0)

Prep PTQ S1
Convert to numeric

s1_ptq_num <- foreach (i = 1:ncol(s1_ptq)) %do% {
  data.frame(as.numeric(unlist(s1_ptq[i]))) %>% setNames(names(s1_ptq[i]))
}%>% bind_cols()

Just 4 data points (< .1% of data) missing for remaining, so mean impute these

length(which(is.na(s1_ptq_num[, 1:15])))
## [1] 4
length(which(is.na(s1_ptq_num[, 1:15])))/(dim(s1_ptq_num[, 1:15])[1]*dim(s1_ptq_num[, 1:15])[2]) 
## [1] 0.002185792
s1_ptq_final <- foreach (i = 1:nrow(s1_ptq_num)) %do% {
  if (any(is.na(s1_ptq_num[i, ]))) {
    s1_ptq_num[i, is.na(s1_ptq_num[i, ])] <- mean(unlist(s1_ptq_num[i, ]), na.rm=TRUE)
  }
s1_ptq_num[i, ]
}%>% bind_rows()
hist(rowSums(s1_ptq_final[, 1:15]), breaks=25)

s1_ptq_final$ptq_sum <- rowSums(s1_ptq_final[, 1:15])
s1_ptq_final$ID <- s1_qs$ID

Sanity check that correlated in expected irection

ComparePars(s1_ptq_final$ptq_sum, s1_bis_final$bis_sum, use_identity_line = 0)

ComparePars(s1_ptq_final$ptq_sum, s1_bas_final$bas_sum, use_identity_line = 0)

Spot check IDs lined up properly

# s1_ptq_final %>% filter(ID == 22)
# s1_ptq_copy %>% filter(ID == 22)
# s1_qs %>% filter(ID == 22)
# s1_ptq_final %>% filter(ID == 104)
# s1_ptq_copy %>% filter(ID == 104)
# s1_qs %>% filter(ID == 104)

Prep S1 RRS

s1_rrs_num <- foreach (i = 1:ncol(s1_rrs)) %do% {
  data.frame(as.numeric(unlist(s1_rrs[i]))) %>% setNames(names(s1_rrs[i]))
}%>% bind_cols()
length(which(is.na(s1_rrs_num[, 1:10])))
## [1] 3
length(which(is.na(s1_rrs_num[, 1:10])))/(dim(s1_rrs_num[, 1:10])[1]*dim(s1_rrs_num[, 1:10])[2]) 
## [1] 0.002459016
s1_rrs_final <- foreach (i = 1:nrow(s1_rrs_num)) %do% {
  if (any(is.na(s1_rrs_num[i, ]))) {
    s1_rrs_num[i, is.na(s1_rrs_num[i, ])] <- mean(unlist(s1_rrs_num[i, ]), na.rm=TRUE)
  }
s1_rrs_num[i, ]
}%>% bind_rows()
hist(rowSums(s1_rrs_final[, 1:10]), breaks=25)

s1_rrs_final$rrs_sum <- rowSums(s1_rrs_final[, 1:10])
s1_rrs_final$ID <- s1_qs$ID

s1 Spot check

# s1_rrs_final %>% filter(ID == 106)
# s1_rrs_copy %>% filter(ID == 106)
# s1_qs %>% filter(ID == 106) # Q8 is RRS

Sanity check RRS and PTQ are strongly correlated

ComparePars(s1_rrs_final$rrs_sum, s1_ptq_final$ptq_sum, use_identity_line = 0)

ComparePars(s1_rrs_final$rrs_sum, s1_bis_final$bis_sum, use_identity_line = 0)

ComparePars(s1_rrs_final$rrs_sum, s1_bas_final$bas_sum, use_identity_line = 0)

Reduce to just matching IDs

m1_s1_model_red <- m1_study1_eb %>% filter(ID %in% s1_ptq_final$ID)
m1_s1_model_red$c_min_o_phi <- m1_s1_model_red$cog_phi - m1_s1_model_red$overt_phi
assert(all(m1_s1_model_red$ID==s1_ptq_final$ID))
assert(all(m1_s1_model_red$ID==s1_rrs_final$ID))

PTQ

ComparePars(m1_s1_model_red$c_min_o_phi, s1_ptq_final$ptq_sum, use_identity_line = 0)

ComparePars(m1_s1_model_red$cog_phi, s1_ptq_final$ptq_sum, use_identity_line = 0)

ComparePars(m1_s1_model_red$overt_phi, s1_ptq_final$ptq_sum, use_identity_line = 0)

Same basic pattern for RRS 4/4 - this was where p value for top statistic slightly different due to imputation error in first version involving ID, adjusted in manuscript from

ComparePars(m1_s1_model_red$c_min_o_phi, s1_rrs_final$rrs_sum, use_identity_line = 0)

# Would expect higher decay to positively correlate w more brooding — instead small and ns neg  
ComparePars(m1_s1_model_red$cog_phi, s1_rrs_final$rrs_sum, use_identity_line = 0)

ComparePars(m1_s1_model_red$overt_phi, s1_rrs_final$rrs_sum, use_identity_line = 0)

Do same thing with age

Complete age data so don’t need to reduce

assert(all(s1_demogs$ID==m1_study1_eb$ID))
m1_study1_eb$c_min_o_phi <- m1_study1_eb$cog_phi - m1_study1_eb$overt_phi

Added for second revision — Age — surprisingly not related to particularly high decay in cog. Closest to signif is actually higher decay in overt!

ComparePars(m1_study1_eb$c_min_o_phi, s1_demogs$Age, use_identity_line = 0)

ComparePars(m1_study1_eb$cog_phi, s1_demogs$Age, use_identity_line = 0)

ComparePars(m1_study1_eb$overt_phi, s1_demogs$Age, use_identity_line = 0)

Age and task performance

s1_train_perf <- s1_train %>% group_by(ID) %>% summarize(m=mean(correct))
s1_test_perf <- s1_sit %>% group_by(ID) %>% summarize(m=mean(correct))
assert(all(s1_demogs$ID==s1_train_perf$ID))
assert(all(s1_demogs$ID==s1_test_perf$ID))

Related to overall performance

ComparePars(s1_train_perf$m, s1_demogs$Age, use_identity_line = 0)

ComparePars(s1_test_perf$m, s1_demogs$Age, use_identity_line = 0)

But not difference in cog vs. overt performance

s1_perf <- s1_train %>% group_by(condition, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
o_min_c_perf_s1 <- as.numeric(unlist(s1_perf[s1_perf$condition=="overt","m"]-
                    s1_perf[s1_perf$condition=="cognitive", "m"]))

# .. and test subject-level differences in performance  
s1_perf_test <- s1_sit %>% group_by(condition, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
o_min_c_testperf_s1 <- as.numeric(unlist(s1_perf_test[s1_perf_test$condition=="overt","m"]-
                    s1_perf_test[s1_perf_test$condition=="cognitive", "m"]))
ComparePars(o_min_c_perf_s1, s1_demogs$Age, use_identity_line = 0)

ComparePars(o_min_c_testperf_s1, s1_demogs$Age, use_identity_line = 0)

Added for second revision — any relationships with reward vs. punishment learning and BIS/BAS?

s1_rew_train_perf <- s1_train %>% filter(valence=="reward") %>% group_by(ID) %>% summarize(m=mean(correct))
s1_pun_train_perf <- s1_train %>% filter(valence=="punishment") %>% group_by(ID) %>% summarize(m=mean(correct))
s1_rew_train_perf$rew_min_pun <-  s1_rew_train_perf$m-s1_pun_train_perf$m

s1_rew_sit_perf <- s1_sit %>% filter(valence=="reward") %>% group_by(ID) %>% summarize(m=mean(correct))
s1_pun_sit_perf <- s1_sit %>% filter(valence=="punishment") %>% group_by(ID) %>% summarize(m=mean(correct))
s1_rew_sit_perf$rew_min_pun <- s1_rew_sit_perf$m-s1_pun_sit_perf$m

ComparePars(s1_rew_train_perf$rew_min_pun, s1_rew_sit_perf$rew_min_pun, use_identity_line = 0)

Reduce to just matching IDs

s1_rew_train_perf_red <- s1_rew_train_perf %>% filter(ID %in% s1_bis_final$ID)
s1_rew_sit_perf_red <- s1_rew_sit_perf %>% filter(ID %in% s1_bis_final$ID)
ComparePars(s1_rew_train_perf_red$rew_min_pun, s1_rew_sit_perf_red$rew_min_pun, use_identity_line = 0)

assert(all(s1_rew_sit_perf_red$ID==s1_bis_final$ID))
assert(all(s1_rew_sit_perf_red$ID==s1_bas_final$ID))
assert(all(s1_rew_train_perf_red$ID==s1_bis_final$ID))
assert(all(s1_rew_train_perf_red$ID==s1_bas_final$ID))

Not correlated

ComparePars(s1_rew_train_perf_red$rew_min_pun, s1_bis_final$bis_sum, use_identity_line = 0)

ComparePars(s1_rew_train_perf_red$rew_min_pun, s1_bas_final$bas_sum, use_identity_line = 0)

ComparePars(s1_rew_sit_perf_red$rew_min_pun, s1_bis_final$bis_sum, use_identity_line = 0)

ComparePars(s1_rew_sit_perf_red$rew_min_pun, s1_bas_final$bas_sum, use_identity_line = 0)

s2 questionnaires new code

No pts with all skips here and table above shows no exclusions

Get and recode PTQ and RRS

s2_ptq <- s2_qs %>% select(contains("Q2_")) %>% setNames(paste0("PTQ_", 1:15))
#s2_ptq$ID <- s2_qs$ID 4/4 removing to make sure not in imputation 
s2_ptq_copy <- s2_ptq
s2_ptq_copy$ID <- s2_qs$ID
s2_ptq[s2_ptq=="Never"] <- 0
s2_ptq[s2_ptq=="Rarely"] <- 1
s2_ptq[s2_ptq=="Sometimes"] <- 2 
s2_ptq[s2_ptq=="Often"] <- 3
s2_ptq[s2_ptq=="Almost always"] <- 4
s2_rrs <- 
  s2_qs %>% select(contains("Q8_")) %>% setNames(paste0("RRS_", 1:10))
#s2_rrs$ID <- s2_qs$ID 4/4 removing to make sure not in imputation 
s2_rrs_copy <- s2_rrs
s2_rrs_copy$ID <- s2_qs$ID
s2_rrs[s2_rrs=="Almost never"] <- 1
s2_rrs[s2_rrs=="Sometimes"] <- 2
s2_rrs[s2_rrs=="Often"] <- 3
s2_rrs[s2_rrs=="Almost always"] <- 4

For revision: BIS/BAS

Exclude catch but keep filler items for now

#s2_qs %>% select(contains("Q7_")) %>% select("Q7_19") # confirmed this is catch 
s2_bis_bas <- s2_qs %>% select(contains("Q7_")) %>% select(-"Q7_19") 

“Items other than 2 and 22 are reverse coded”

The two BIS items that are forward coded (1= very true, 4=very false) — Q7_23 corresponds to item 22 in actual scale because of the catch item

fwd_code_bis_bas <- s2_bis_bas %>% select("Q7_2", "Q7_23")
fwd_code_bis_bas[fwd_code_bis_bas=="Very true for me"] <- 1
fwd_code_bis_bas[fwd_code_bis_bas=="Somewhat true for me"] <- 2
fwd_code_bis_bas[fwd_code_bis_bas=="Somewhat false for me"] <- 3
fwd_code_bis_bas[fwd_code_bis_bas=="Very false for me"] <- 4

Sanity check that positively correlate

cor.test(as.numeric(fwd_code_bis_bas$Q7_2), as.numeric(fwd_code_bis_bas$Q7_23))
## 
##  Pearson's product-moment correlation
## 
## data:  as.numeric(fwd_code_bis_bas$Q7_2) and as.numeric(fwd_code_bis_bas$Q7_23)
## t = 8.41, df = 134, p-value = 5.365e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4654952 0.6880297
## sample estimates:
##       cor 
## 0.5877714

The rest are reverse coded

rev_code_bis_bas <- s2_bis_bas %>% select(-c("Q7_2", "Q7_23"))
rev_code_bis_bas[rev_code_bis_bas=="Very true for me"] <- 4
rev_code_bis_bas[rev_code_bis_bas=="Somewhat true for me"] <- 3
rev_code_bis_bas[rev_code_bis_bas=="Somewhat false for me"] <- 2
rev_code_bis_bas[rev_code_bis_bas=="Very false for me"] <- 1
s2_bis_bas <- data.frame(fwd_code_bis_bas, rev_code_bis_bas)
s2_bis_bas_num <- foreach (i = 1:ncol(s2_bis_bas)) %do% {
  data.frame(as.numeric(unlist(s2_bis_bas[i]))) %>% setNames(names(s2_bis_bas[i]))
}%>% bind_cols()

Sanity check BIS fwd and reverse code items positively correlate

cor.test(s2_bis_bas_num$Q7_2, s2_bis_bas_num$Q7_23)
## 
##  Pearson's product-moment correlation
## 
## data:  s2_bis_bas_num$Q7_2 and s2_bis_bas_num$Q7_23
## t = 8.41, df = 134, p-value = 5.365e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4654952 0.6880297
## sample estimates:
##       cor 
## 0.5877714
cor.test(s2_bis_bas_num$Q7_2, s2_bis_bas_num$Q7_8)
## 
##  Pearson's product-moment correlation
## 
## data:  s2_bis_bas_num$Q7_2 and s2_bis_bas_num$Q7_8
## t = 9.1229, df = 134, p-value = 9.631e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5030616 0.7130186
## sample estimates:
##       cor 
## 0.6189783
cor.test(s2_bis_bas_num$Q7_23, s2_bis_bas_num$Q7_8)
## 
##  Pearson's product-moment correlation
## 
## data:  s2_bis_bas_num$Q7_23 and s2_bis_bas_num$Q7_8
## t = 6.2383, df = 135, p-value = 5.312e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3316317 0.5936545
## sample estimates:
##      cor 
## 0.473037
# Item 24 because shifted by 1 due to catch
cor.test(s2_bis_bas_num$Q7_2, s2_bis_bas_num$Q7_25)
## 
##  Pearson's product-moment correlation
## 
## data:  s2_bis_bas_num$Q7_2 and s2_bis_bas_num$Q7_25
## t = 10.117, df = 134, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5507643 0.7439995
## sample estimates:
##       cor 
## 0.6580849
cor.test(s2_bis_bas_num$Q7_23, s2_bis_bas_num$Q7_25)
## 
##  Pearson's product-moment correlation
## 
## data:  s2_bis_bas_num$Q7_23 and s2_bis_bas_num$Q7_25
## t = 6.5828, df = 135, p-value = 9.4e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3545357 0.6102077
## sample estimates:
##       cor 
## 0.4929404
cor.test(s2_bis_bas_num$Q7_2, s2_bis_bas_num$Q7_16)
## 
##  Pearson's product-moment correlation
## 
## data:  s2_bis_bas_num$Q7_2 and s2_bis_bas_num$Q7_16
## t = 12.059, df = 134, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6295342 0.7934011
## sample estimates:
##       cor 
## 0.7214175
cor.test(s2_bis_bas_num$Q7_23, s2_bis_bas_num$Q7_16)
## 
##  Pearson's product-moment correlation
## 
## data:  s2_bis_bas_num$Q7_23 and s2_bis_bas_num$Q7_16
## t = 8.874, df = 135, p-value = 3.773e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4890406 0.7031120
## sample estimates:
##       cor 
## 0.6069724

Now remove filler items — all before catch so all items corresopnd to actual num in scale

s2_bis_bas <- s2_bis_bas_num %>% select(-c(
  "Q7_1",
  "Q7_6",
  "Q7_11",
  "Q7_17") 
)

BIS

s2_bis <- s2_bis_bas %>% select(
  # Before catch item 
  c("Q7_2", "Q7_8", "Q7_13", "Q7_16"),
  # After catch so one shifted
  c("Q7_20", "Q7_23", "Q7_25"))

Just 1 missing so impute

length(which(is.na(s2_bis[, 1:7])))
## [1] 1
length(which(is.na(s2_bis[, 1:7])))/(dim(s2_bis[, 1:7])[1]*dim(s2_bis[, 1:7])[2]) 
## [1] 0.001042753
s2_bis_final <- foreach (i = 1:nrow(s2_bis)) %do% {
  if (any(is.na(s2_bis[i, ]))) {
    s2_bis[i, is.na(s2_bis[i, ])] <- mean(unlist(s2_bis[i, ]), na.rm=TRUE)
  }
s2_bis[i, ]
}%>% bind_rows()

Sanity check all positively correlated

cor(s2_bis_final)
##            Q7_2      Q7_8     Q7_13     Q7_16     Q7_20     Q7_23     Q7_25
## Q7_2  1.0000000 0.6189294 0.5450794 0.7215200 0.5609502 0.5875614 0.6590967
## Q7_8  0.6189294 1.0000000 0.6954469 0.6496518 0.6489944 0.4730370 0.7506610
## Q7_13 0.5450794 0.6954469 1.0000000 0.6336202 0.7315979 0.3881398 0.6686291
## Q7_16 0.7215200 0.6496518 0.6336202 1.0000000 0.6550656 0.6069724 0.7005364
## Q7_20 0.5609502 0.6489944 0.7315979 0.6550656 1.0000000 0.3751194 0.7093946
## Q7_23 0.5875614 0.4730370 0.3881398 0.6069724 0.3751194 1.0000000 0.4929404
## Q7_25 0.6590967 0.7506610 0.6686291 0.7005364 0.7093946 0.4929404 1.0000000

And high alpha

psych::alpha(s2_bis_final)
## 
## Reliability analysis   
## Call: psych::alpha(x = s2_bis_final)
## 
##   raw_alpha std.alpha G6(smc) average_r S/N   ase mean   sd median_r
##       0.92      0.92    0.92      0.61  11 0.011    3 0.76     0.65
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt     0.89  0.92  0.94
## Duhachek  0.90  0.92  0.94
## 
##  Reliability if an item is dropped:
##       raw_alpha std.alpha G6(smc) average_r  S/N alpha se  var.r med.r
## Q7_2       0.91      0.90    0.90      0.61  9.5    0.012 0.0147  0.65
## Q7_8       0.90      0.90    0.90      0.60  9.1    0.013 0.0127  0.63
## Q7_13      0.91      0.91    0.90      0.61  9.5    0.012 0.0106  0.65
## Q7_16      0.90      0.90    0.90      0.59  8.8    0.013 0.0142  0.62
## Q7_20      0.91      0.90    0.90      0.61  9.5    0.012 0.0101  0.63
## Q7_23      0.92      0.92    0.92      0.66 11.8    0.010 0.0034  0.66
## Q7_25      0.90      0.90    0.90      0.59  8.7    0.013 0.0119  0.62
## 
##  Item statistics 
##         n raw.r std.r r.cor r.drop mean   sd
## Q7_2  137  0.82  0.82  0.78   0.75  3.1 0.91
## Q7_8  137  0.85  0.85  0.82   0.78  2.8 0.99
## Q7_13 137  0.82  0.81  0.78   0.74  3.0 0.99
## Q7_16 137  0.87  0.87  0.85   0.81  2.7 0.98
## Q7_20 137  0.82  0.82  0.79   0.75  3.2 0.86
## Q7_23 137  0.68  0.69  0.61   0.58  2.9 0.85
## Q7_25 137  0.87  0.87  0.85   0.82  3.0 0.95
## 
## Non missing response frequency for each item
##          1    2    3 3.66666666666667    4 miss
## Q7_2  0.06 0.18 0.33             0.01 0.43    0
## Q7_8  0.11 0.26 0.33             0.00 0.30    0
## Q7_13 0.09 0.21 0.31             0.00 0.39    0
## Q7_16 0.13 0.26 0.37             0.00 0.24    0
## Q7_20 0.05 0.12 0.38             0.00 0.45    0
## Q7_23 0.04 0.29 0.39             0.00 0.27    0
## Q7_25 0.09 0.17 0.36             0.00 0.39    0
s2_bis_final$bis_sum <- rowSums(s2_bis_final)
s2_bis_final$ID <- s2_qs$ID

The remaining items are BAS

s2_bas <- s2_bis_bas %>% select(
  
  -c("Q7_2", "Q7_8", "Q7_13", "Q7_16",
  
  "Q7_20", "Q7_23", "Q7_25"))

Nothing missing

s2_bas_final <- s2_bas

Sanity check all positively correlated

cor(s2_bas_final)
##            Q7_3      Q7_4      Q7_5      Q7_7      Q7_9     Q7_10     Q7_12
## Q7_3  1.0000000 0.3185535 0.3257158 0.4426593 0.5869646 0.1445450 0.5553761
## Q7_4  0.3185535 1.0000000 0.3380587 0.5274478 0.2944064 0.1264353 0.3069515
## Q7_5  0.3257158 0.3380587 1.0000000 0.4129971 0.3148678 0.4485842 0.2683164
## Q7_7  0.4426593 0.5274478 0.4129971 1.0000000 0.3433022 0.2334689 0.3594777
## Q7_9  0.5869646 0.2944064 0.3148678 0.3433022 1.0000000 0.2683965 0.6989141
## Q7_10 0.1445450 0.1264353 0.4485842 0.2334689 0.2683965 1.0000000 0.2664901
## Q7_12 0.5553761 0.3069515 0.2683164 0.3594777 0.6989141 0.2664901 1.0000000
## Q7_14 0.3922082 0.4489550 0.4300696 0.5804649 0.5256589 0.2812001 0.5457326
## Q7_15 0.1813896 0.1093602 0.2558661 0.2474474 0.3159382 0.5161126 0.2604770
## Q7_18 0.2252178 0.3198891 0.2027570 0.3412878 0.2780012 0.2329042 0.2405645
## Q7_21 0.2577641 0.2007712 0.5532761 0.2822111 0.2790424 0.3197975 0.2602531
## Q7_22 0.3994105 0.1852181 0.2325790 0.1978012 0.5694647 0.2151834 0.4634925
## Q7_24 0.2786192 0.3872116 0.3449801 0.5390740 0.2116188 0.2523446 0.2028212
##           Q7_14      Q7_15      Q7_18     Q7_21     Q7_22     Q7_24
## Q7_3  0.3922082 0.18138960 0.22521779 0.2577641 0.3994105 0.2786192
## Q7_4  0.4489550 0.10936015 0.31988905 0.2007712 0.1852181 0.3872116
## Q7_5  0.4300696 0.25586613 0.20275702 0.5532761 0.2325790 0.3449801
## Q7_7  0.5804649 0.24744745 0.34128782 0.2822111 0.1978012 0.5390740
## Q7_9  0.5256589 0.31593825 0.27800116 0.2790424 0.5694647 0.2116188
## Q7_10 0.2812001 0.51611265 0.23290420 0.3197975 0.2151834 0.2523446
## Q7_12 0.5457326 0.26047702 0.24056449 0.2602531 0.4634925 0.2028212
## Q7_14 1.0000000 0.21919045 0.42229892 0.3227852 0.2884935 0.4626583
## Q7_15 0.2191904 1.00000000 0.09816021 0.3793122 0.3428206 0.1690617
## Q7_18 0.4222989 0.09816021 1.00000000 0.2490459 0.3142644 0.3727922
## Q7_21 0.3227852 0.37931218 0.24904588 1.0000000 0.3303116 0.2339747
## Q7_22 0.2884935 0.34282056 0.31426445 0.3303116 1.0000000 0.1095123
## Q7_24 0.4626583 0.16906170 0.37279220 0.2339747 0.1095123 1.0000000

And high alpha

psych::alpha(s2_bas_final) 
## 
## Reliability analysis   
## Call: psych::alpha(x = s2_bas_final)
## 
##   raw_alpha std.alpha G6(smc) average_r S/N   ase mean   sd median_r
##       0.86      0.86    0.89      0.33 6.3 0.017  2.9 0.46     0.31
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt     0.82  0.86  0.89
## Duhachek  0.83  0.86  0.90
## 
##  Reliability if an item is dropped:
##       raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## Q7_3       0.85      0.85    0.88      0.32 5.7    0.019 0.016  0.30
## Q7_4       0.86      0.86    0.88      0.33 6.0    0.018 0.017  0.30
## Q7_5       0.85      0.85    0.87      0.32 5.7    0.019 0.018  0.29
## Q7_7       0.85      0.85    0.87      0.32 5.6    0.019 0.016  0.29
## Q7_9       0.84      0.85    0.87      0.31 5.5    0.020 0.014  0.30
## Q7_10      0.86      0.86    0.88      0.34 6.1    0.018 0.017  0.32
## Q7_12      0.85      0.85    0.87      0.32 5.6    0.019 0.015  0.31
## Q7_14      0.84      0.84    0.87      0.31 5.4    0.020 0.016  0.28
## Q7_15      0.86      0.86    0.88      0.34 6.2    0.018 0.016  0.32
## Q7_18      0.86      0.86    0.88      0.34 6.1    0.018 0.018  0.32
## Q7_21      0.85      0.86    0.88      0.33 5.9    0.018 0.018  0.31
## Q7_22      0.85      0.86    0.88      0.33 5.9    0.018 0.017  0.31
## Q7_24      0.86      0.86    0.88      0.33 6.0    0.018 0.017  0.31
## 
##  Item statistics 
##         n raw.r std.r r.cor r.drop mean   sd
## Q7_3  137  0.64  0.64  0.61   0.56  2.7 0.76
## Q7_4  137  0.54  0.57  0.52   0.47  3.5 0.61
## Q7_5  137  0.65  0.64  0.61   0.55  3.0 0.82
## Q7_7  137  0.67  0.69  0.67   0.59  3.4 0.70
## Q7_9  137  0.73  0.71  0.70   0.65  2.7 0.84
## Q7_10 137  0.55  0.54  0.49   0.45  2.7 0.81
## Q7_12 137  0.69  0.68  0.66   0.61  2.7 0.79
## Q7_14 137  0.73  0.74  0.73   0.67  3.1 0.73
## Q7_15 137  0.52  0.51  0.46   0.42  2.1 0.76
## Q7_18 137  0.53  0.54  0.48   0.43  3.0 0.80
## Q7_21 137  0.59  0.58  0.54   0.50  2.6 0.79
## Q7_22 137  0.60  0.58  0.54   0.50  2.2 0.83
## Q7_24 137  0.54  0.57  0.52   0.47  3.5 0.57
## 
## Non missing response frequency for each item
##          1    2    3    4 miss
## Q7_3  0.05 0.31 0.50 0.13    0
## Q7_4  0.00 0.06 0.37 0.57    0
## Q7_5  0.04 0.20 0.46 0.30    0
## Q7_7  0.01 0.10 0.39 0.50    0
## Q7_9  0.08 0.34 0.42 0.16    0
## Q7_10 0.08 0.28 0.50 0.14    0
## Q7_12 0.06 0.31 0.48 0.15    0
## Q7_14 0.02 0.17 0.54 0.27    0
## Q7_15 0.23 0.53 0.21 0.04    0
## Q7_18 0.04 0.21 0.48 0.27    0
## Q7_21 0.07 0.38 0.43 0.12    0
## Q7_22 0.22 0.46 0.26 0.06    0
## Q7_24 0.00 0.04 0.43 0.53    0
s2_bas_final$bas_sum <- rowSums(s2_bas_final)
hist(s2_bas_final$bas_sum, breaks=25)

s2_bas_final$ID <- s2_qs$ID

Modestly inversely related

ComparePars(s2_bas_final$bas_sum, s2_bis_final$bis_sum, use_identity_line = 0)

Prep PTQ s2
Convert to numeric

s2_ptq_num <- foreach (i = 1:ncol(s2_ptq)) %do% {
  data.frame(as.numeric(unlist(s2_ptq[i]))) %>% setNames(names(s2_ptq[i]))
}%>% bind_cols()

Just 3 data points (< .1% of data) missing for remaining, so mean impute these

length(which(is.na(s2_ptq_num[, 1:15])))
## [1] 3
length(which(is.na(s2_ptq_num[, 1:15])))/(dim(s2_ptq_num[, 1:15])[1]*dim(s2_ptq_num[, 1:15])[2]) 
## [1] 0.001459854
s2_ptq_final <- foreach (i = 1:nrow(s2_ptq_num)) %do% {
  if (any(is.na(s2_ptq_num[i, ]))) {
    s2_ptq_num[i, is.na(s2_ptq_num[i, ])] <- mean(unlist(s2_ptq_num[i, ]), na.rm=TRUE)
  }
s2_ptq_num[i, ]
}%>% bind_rows()
hist(rowSums(s2_ptq_final[, 1:15]), breaks=25)

s2_ptq_final$ptq_sum <- rowSums(s2_ptq_final[, 1:15])
s2_ptq_final$ID <- s2_qs$ID

Sanity check that correlated in expected direction

ComparePars(s2_ptq_final$ptq_sum, s2_bis_final$bis_sum, use_identity_line = 0)

ComparePars(s2_ptq_final$ptq_sum, s2_bas_final$bas_sum, use_identity_line = 0)

Spot check IDs lined up properly

# s2_ptq_final %>% filter(ID == 22)
# s2_ptq_copy %>% filter(ID == 22)
# # s2_qs %>% filter(ID == 22)
# s2_ptq_final %>% filter(ID == 104)
# s2_ptq_copy %>% filter(ID == 104)
# s2_qs %>% filter(ID == 104)

Prep s2 RRS

s2_rrs_num <- foreach (i = 1:ncol(s2_rrs)) %do% {
  data.frame(as.numeric(unlist(s2_rrs[i]))) %>% setNames(names(s2_rrs[i]))
}%>% bind_cols()
length(which(is.na(s2_rrs_num[, 1:10])))
## [1] 3
length(which(is.na(s2_rrs_num[, 1:10])))/(dim(s2_rrs_num[, 1:10])[1]*dim(s2_rrs_num[, 1:10])[2]) 
## [1] 0.002189781
s2_rrs_final <- foreach (i = 1:nrow(s2_rrs_num)) %do% {
  if (any(is.na(s2_rrs_num[i, ]))) {
    s2_rrs_num[i, is.na(s2_rrs_num[i, ])] <- mean(unlist(s2_rrs_num[i, ]), na.rm=TRUE)
  }
s2_rrs_num[i, ]
}%>% bind_rows()
hist(rowSums(s2_rrs_final[, 1:10]), breaks=25)

s2_rrs_final$rrs_sum <- rowSums(s2_rrs_final[, 1:10])
s2_rrs_final$ID <- s2_qs$ID

s2 Spot check

# s2_rrs_final %>% filter(ID == 106)
# s2_rrs_copy %>% filter(ID == 106)
# s2_qs %>% filter(ID == 106) # Q8 is RRS

Sanity check RRS and PTQ are strongly correlated

ComparePars(s2_rrs_final$rrs_sum, s2_ptq_final$ptq_sum, use_identity_line = 0)

ComparePars(s2_rrs_final$rrs_sum, s2_bis_final$bis_sum, use_identity_line = 0)

ComparePars(s2_rrs_final$rrs_sum, s2_bas_final$bas_sum, use_identity_line = 0)

Reduce to just matching IDs

m1_s2_model_red <- m1_study2_eb %>% filter(ID %in% s2_ptq_final$ID)
m1_s2_model_red$c_min_o_phi <- m1_s2_model_red$cog_phi - m1_s2_model_red$overt_phi
assert(all(m1_s2_model_red$ID==s2_ptq_final$ID))
assert(all(m1_s2_model_red$ID==s2_rrs_final$ID))

PTQ

ComparePars(m1_s2_model_red$c_min_o_phi, s2_ptq_final$ptq_sum, use_identity_line = 0)

ComparePars(m1_s2_model_red$cog_phi, s2_ptq_final$ptq_sum, use_identity_line = 0)

ComparePars(m1_s2_model_red$overt_phi, s2_ptq_final$ptq_sum, use_identity_line = 0)

Same basic pattern for RRS

ComparePars(m1_s2_model_red$c_min_o_phi, s2_rrs_final$rrs_sum, use_identity_line = 0)

# Would expect higher decay to positively correlate w more brooding — instead small and ns neg  
ComparePars(m1_s2_model_red$cog_phi, s2_rrs_final$rrs_sum, use_identity_line = 0)

ComparePars(m1_s2_model_red$overt_phi, s2_rrs_final$rrs_sum, use_identity_line = 0)

s2_rew_train_perf <- s2_train %>% filter(valence=="reward") %>% group_by(ID) %>% summarize(m=mean(correct))
s2_pun_train_perf <- s2_train %>% filter(valence=="punishment") %>% group_by(ID) %>% summarize(m=mean(correct))
s2_rew_train_perf$rew_min_pun <-  s2_rew_train_perf$m-s2_pun_train_perf$m

s2_rew_sit_perf <- s2_sit %>% filter(valence=="reward") %>% group_by(ID) %>% summarize(m=mean(correct))
s2_pun_sit_perf <- s2_sit %>% filter(valence=="punishment") %>% group_by(ID) %>% summarize(m=mean(correct))
s2_rew_sit_perf$rew_min_pun <- s2_rew_sit_perf$m-s2_pun_sit_perf$m

ComparePars(s2_rew_train_perf$rew_min_pun, s2_rew_sit_perf$rew_min_pun, use_identity_line = 0)

Reduce to just matching IDs

s2_rew_train_perf_red <- s2_rew_train_perf %>% filter(ID %in% s2_bis_final$ID)
s2_rew_sit_perf_red <- s2_rew_sit_perf %>% filter(ID %in% s2_bis_final$ID)
ComparePars(s2_rew_train_perf_red$rew_min_pun, s2_rew_sit_perf_red$rew_min_pun, use_identity_line = 0)

assert(all(s2_rew_sit_perf_red$ID==s2_bis_final$ID))
assert(all(s2_rew_sit_perf_red$ID==s2_bas_final$ID))
assert(all(s2_rew_train_perf_red$ID==s2_bis_final$ID))
assert(all(s2_rew_train_perf_red$ID==s2_bas_final$ID))

Not correlated

ComparePars(s2_rew_train_perf_red$rew_min_pun, s2_bis_final$bis_sum, use_identity_line = 0)

ComparePars(s2_rew_train_perf_red$rew_min_pun, s2_bas_final$bas_sum, use_identity_line = 0)

ComparePars(s2_rew_sit_perf_red$rew_min_pun, s2_bis_final$bis_sum, use_identity_line = 0)

ComparePars(s2_rew_sit_perf_red$rew_min_pun, s2_bas_final$bas_sum, use_identity_line = 0)

Added for second revision — Age

There are some Prolific pts without this so take them out and match the datasets

s2_demogs_red <- s2_demogs %>% filter(Age != "DATA_EXPIRED")
s2_demogs_red$Age <- as.numeric(s2_demogs_red$Age)
m1_study2_age_red <- m1_study2_eb %>% filter(ID %in% s2_demogs_red$ID)

assert(all(s2_demogs_red$ID==m1_study2_age_red$ID))
m1_study2_age_red$c_min_o_phi <- m1_study2_age_red$cog_phi - m1_study2_age_red$overt_phi

Only one significant is overt decay!

ComparePars(m1_study2_age_red$c_min_o_phi, s2_demogs_red$Age, use_identity_line = 0)

ComparePars(m1_study2_age_red$cog_phi, s2_demogs_red$Age, use_identity_line = 0)

ComparePars(m1_study2_age_red$overt_phi, s2_demogs_red$Age, use_identity_line = 0)

Age and task performance

s2_train_perf_age_red <- s2_train %>% filter(ID %in% s2_demogs_red$ID) %>% group_by(ID) %>% summarize(m=mean(correct))
s2_test_perf_age_red <- s2_sit %>% filter(ID %in% s2_demogs_red$ID) %>%  group_by(ID) %>% summarize(m=mean(correct))
assert(all(s2_demogs_red$ID==s2_train_perf_age_red$ID))
assert(all(s2_demogs_red$ID==s2_test_perf_age_red$ID))

In study 2 unrelated to overall performance

ComparePars(s2_train_perf_age_red$m, s2_demogs_red$Age, use_identity_line = 0)

ComparePars(s2_test_perf_age_red$m, s2_demogs_red$Age, use_identity_line = 0)

s2_perf_age_red <- s2_train %>% filter(ID %in% s2_demogs_red$ID) %>%
  group_by(condition, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
o_min_c_perf_s2 <- as.numeric(unlist(s2_perf_age_red[s2_perf_age_red$condition=="overt","m"]-
                    s2_perf_age_red[s2_perf_age_red$condition=="cognitive", "m"]))

# .. and test subject-level differences in performance  
s2_perf_age_red_test <- s2_sit %>% filter(ID %in% s2_demogs_red$ID) %>% 
  group_by(condition, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
o_min_c_testperf_s2 <- as.numeric(unlist(s2_perf_age_red_test[s2_perf_age_red_test$condition=="overt","m"]-
                    s2_perf_age_red_test[s2_perf_age_red_test$condition=="cognitive", "m"]))
ComparePars(o_min_c_perf_s2, s2_demogs_red$Age, use_identity_line = 0)

ComparePars(o_min_c_testperf_s2, s2_demogs_red$Age, use_identity_line = 0)